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ABSTRACT 


In an effort to pursue advanced space missions, improved onboard trajectory opti¬ 
mization and path-planning capabilities are necessary. Regardless of the mission, the 
paramount requirement for any candidate guidance-and-control method is the ability 
to react in real time to a dynamic environment, followed by fuel efficiency. Ground- 
based kinodynamic test beds are critical in developing, testing, and verifying these 
requirements. The NFS POSEIDYN test bed is introduced here as a state-of-the-art 
(SOTA) dynamic, hardware-in-the-loop test bed. Key improvements to the software 
architecture, enabling the development of multi-rate guidance, navigation, and control 
(GNC) algorithms, are presented in addition to a detailed system characterization. To 
aid in the experimental evaluation of GNC algorithms, a Standard Test Framework is 
proposed in addition to a new guidance comparison metric previously missing in the 
literature, which simultaneously captures the computational complexity and algorithm 
performance. Guidance and control methods representing the SOTA are experimen¬ 
tally evaluated and compared against a proposed rapidly exploring random tree-based 
guidance method using the proposed test framework and the comparison metric. 
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CHAPTER 1: 
Introduction 


The formulation of the problem is often more essential than its 
solution, which may be merely a matter of mathematical or 
experimental skill. To raise new questions, new possibilities, to 
regard old problems from a new angle requires creative 
imagination and marks real advance in science. 

—Albert Einstein and Leopold Infield 
The Evolution of Physics [1 ] 


1.1 Motivation 

Suppose for a moment, you are tasked with designing the guidance and control (G&C) 
architecture for an interplanetary sample-return mission. For this notional mission, the 
spacecraft will rendezvous with an object, travelling through the expanse of the Solar 
System, that is believed to be (partially) composed of some rare and highly-desired 
material. After it has rendezvoused, it will perform proximity maneuvers to collect vari¬ 
ous measurements as well as a physical sample from the object. After completing the 
mission, the spacecraft will rendezvous and dock with a pre-positioned fuel depot. The 
question becomes, what guidance and control methods will you choose? Why? 

Do the intrinsic characteristics and properties of your selected method make it an at¬ 
tractive choice? Were there any environmental factors - such as communication de¬ 
lays and other naturally occurring bodies (e.g., micrometeoroids) - that require an au¬ 
tonomous response from the spacecraft? If these are important factors, is the original 
G&C still the best choice? Moreover, after you select a method, how are you going to 
test and validate the implementation? 

While this is just a notional thought experiment, it illustrates some of the needs pre¬ 
sented in the National Research Council’s decadal study titled Voyages for Planetary 
Science in the Decade 2013-2022 [2] and technologies identified in the National Aero¬ 
nautics and Space Administration (NASA) Technology Roadmaps [3]. Both studies 
emphasize the need for improvements to onboard autonomous trajectory optimization, 
path planning, and replanning capabilities to support envisioned space missions [2]-[4]. 
These necessary advancements in onboard autonomous Guidance, Navigation, and 
Control (GNC) capabilities will enable more complex missions, such as autonomous 
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spacecraft rendezvous and proximity operations (RPO) as well as autonomous ren¬ 
dezvous and docking. 


1995 


2015 




July 1998 

ETS-VII 

First autonomous 
rendezvous and docking 
between two unmanned 
spacecraft 


Jan 2003 

XSS-10 

Autonomous navigation, , 
prox-ops, and inspection of 
another space object 


March 2008 

ATV 

Autonomous rendezvous 
and docking to resupply 
the ISS 



June 2010 

PRISMA 


Autonomous RPO and tech ^ 
demo of RPO sensors and 
actuators 


April 2005 

XSS-11 


Rendezvous and prox- 
ops and autonomous 
mission planning 




1968 - Present 

Soyuz 

Radio-based navigation 
system (Igla/Kursj used for 
rendezvous and docking 


April 2005 

DART 

Testing of the AVGS 
for relative position 
and orientation 



March 2007 

Orbital Express 

Fully autonomous on-orbit 
rendezvous and servicing 



Dec 2015 - Present 

Falcon 9 
First Stage 

Real-time convex optimization 
for the return of an orbital-class 
reusable launch vehicle stages 



Fall 2017/Spring 2018 

Prox-1 

Fully passive relative navigation 
and first on-orbit use of APF-based 
guidance 



Figure 1.1. Evolution of Real-Time Autonomous Guidance and Control.^ 


As illustrated in Figure 1.1, numerous missions, which were either proposed or flower 
over the past few decades, have explored new techniques for safely conducting RPO 
[5], [7]. Arguably, RPO is considered to have begun in the early 1960s with the 
Vostok^ and Gemini programs while autonomous RPO begin in 1967 with the Soyuz 
program [5]. The more recent push towards greater autonomy to enable more com¬ 
plex missions began with the Japan Aerospace Exploration Agency’s Engineering Test 
Satellite VII (ETS-VII) [5]. Since then, several experimental missions have been com¬ 
pleted (or proposed in the case of the Georgia Institute of Technology’s Prox-1 [7], 
[8]), which have tested various algorithms and equipment necessary to conduct au¬ 
tonomous RPO [5], [9]-[11]. Yet despite active research in this area, few spacecraft 
can be found in literature which routinely perform autonomous RPO and docking. Cur¬ 
rently, the only three spacecraft which autonomously dock with the International Space 
Station (ISS) are the European Space Agency’s Automated Transfer Vehicle (ATV) [18] 
and the Russian Soyuz and Progress spacecraft [5]. Additionally, while not technically 

^The data used to create this figure was adapted from [5H11]; the images were adapted from the 
following sources: [9], [10], [12H17]. 

^The rendezvous was performed by launching a manned capsule in such an orbit that it approached 
another on-orbit capsule (launched the previous day) to within 6.5 km [5]. 
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a spacecraft, the Space Exploration Technologies Corporation (also known as SpaceX) 
has demonstrated the successful use of real-time convex optimization in recovering the 
first stage of its Falcon 9 launch vehicle [6] in an arguably more complex and dynamic 
environment than on orbit — effectively showing the potential for use on-orbit. Regard¬ 
less of the mission, attention must be given to ensure safe operations in the vicinity of 
other spacecraft. 

To achieve autonomous RPO, the guidance algorithm must not only ensure real-time, 
collision-free operations, but also reduce propellant consumption. This results in a 
conflicting set of requirements [19], [20] and is illustrated in Figure 1.2. For instance, 
obstacle avoidance requires the GNC subsystem to have the ability to respond and 
re-plan to a dynamic environment. This is critical, as collision avoidance is a hard 
requirement for any mission [7], [8], [21]. However, generating an efficient and safe 
trajectory involves solving a constrained optimal control problem (OCR). Once the 
optimal solution to the OCR is found, a reference trajectory can then be generated 
and tracked using a series of waypoints [18], [22] or feed-forward control [23]. Given 
the dimensionality of the OCP and possible numerical sensitivities, methods for finding 
optimal trajectories can be too computationally burdensome to be performed onboard 
a computationally constrained spacecraft. Thus, the generation of a sub-optimal but 
feasible solution is more important, as it will ensure the safety-of-fight requirement (i.e., 
collision avoidance) is satisfied at all times. 


Increased Optimalit 
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Figure 1.2. Guidance algorithm trade-off between increased optimality 
and computational burden. 
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1.2 State of the Art Real-Time Guidance and Control 
Methods 

Real-Time Guidance & Control (RT-G&C) algorithms are those that can be implemented 
in the computationally constrained environment present onboard spacecraft and can be 
run in real-time such that they can provide feedback control. By real-time, it is meant 
that the algorithm readily provides a feasible solution within a prescribed amount of 
time regardless of the system load, and failure to meet these deadlines may result 
in a hazardous situation (e.g., a collision with another object) [24], Notably, this time 
period can range from milliseconds to several seconds or minutes and is dependent 
upon several factors (and is further discussed in Section 7.3.1). Real-time systems are 
categorized as soft, hard, and firm [24]. The consequences of missing deadlines as well 
as corresponding examples for each of these categories are presented in Table 1.1. 


Table 1.1. Summary of real-time system categories. Adapted from [24]. 


Category 

Consequence of 

Missed Deadline(s) 

Example Application 

Soft 

Degraded Performance 

Computer Programs 
(e.g.. Word Processor) 

Hard 

Entire System Failure 

Tolerant of the 

Safety-Oriented Systems 
(e.g.. Fire Suppression System) 

Firm 

Occasional Missed Deadline; 
Too Many Missed Deadlines 
May Result in System Failure 

Embedded Controllers 
(e.g.. Engine Management Controller) 


The limiting factor in RT-G&C algorithms is the computational-constrained environment 
present onboard spacecraft. Despite modern computational technology, space-rated 
compute elements are subject to typical design constraints, such as high reliability in 
the presence of radiation, low power consumption, low thermal design power, and high 
computational power [25]. In general, power, volume, mass, and cost constraints make 
adding a dedicated processor for advanced G&C algorithms infeasible. Combined with 
the need to be respect time deadlines, RT-G&C algorithms are inherently limited in their 
computational complexity. 

Different theoretical approaches have been explored to solve the guidance and con¬ 
trol problem for achieving autonomous real-time, collision-free RPO in literature. Note, 
Figure 1.2 lists several state-of-the-art RT-G&C methods that will be discussed. Early 
RT-G&C methods, utilized by re-entry autopilot and reaction-control system in the Apollo 
spacecraft, include phase-plane logic and sliding mode control [26]. More recently. 
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proportional-integral-derivative as well as H-infinity controllers have been utilized by 
the ATV to track a reference trajectory and autonomously dock with the ISS [18]. Other 
approaches in literature are limited to experimental testing only, mostly on terrestrial 
test beds. 

Artificial Potential Function (APF)-based guidance methods have the advantage of 
guaranteeing real-time on-board execution and collision-free RPO through there closed- 
form control law [19], [20], [27], [28]. Additionally, onboard passive or active sensors 
can be used for relative navigation, thus allowing the APF guidance method to achieve 
real-time obstacle avoidance. This method has been successfully implemented for 
autonomous RPO in either simulated or experimental environments [29]-[33]. The per¬ 
formance of the APF method is not guaranteed since there is no explicit minimization 
of a cost function, such as fuel consumption or maneuvering time. 

To improve the performance of APF, a few approaches have been proposed [19], 
[20], [30], [34]. One such approach, the Adaptive Artificial Potential Function (AAPF) 
method, improves the optimality (with respect to a given cost functional) of the APF 
method by adapting the attractive potential shaping matrix to maintain a reference tra¬ 
jectory through the use of an adaptive update law [19], [20]. Note, the adaptive com¬ 
ponent of the AAPF method makes this a hybrid approach, as illustrated by Figure 1.2. 
The AAPF method was shown by Munoz and Fitz-Coy to achieve improved consis¬ 
tency in maneuvering time and reduced propellant usage compared to the APF method 
via Monte Carlo analysis [19], [20]. However, all APF-based methods suffer from the 
existence of (stable) local minima which may trap a spacecraft and prevent it from 
completing the maneuver [20], [26]-[30], [34], [35]. 

To overcome the local minima issue, the use of potential fields based on analytic har¬ 
monic functions was explored by Akishita et al. Harmonic Potential Functions (HPFs) 
use analytical functions that have strong ties to irrotational and inviscid fluid flow; how¬ 
ever, technically speaking, they do not completely eliminate the presence of local min¬ 
ima [36]. Any local minima present in the resulting potential field are unstable (i.e., 
saddle points) [36]-[38]. However, one drawback of HPFs is the difficultly associated 
with properly scaling the potential field elements to reduce numerical noise in the gra¬ 
dient computation [39], [40]. Despite this difficultly, HPFs can be utilized to generate 
a “Navigation Function,” as defined by Rimon and Koditschek [41], due to their con¬ 
struction and analytical formulation. HPF-based methods have been primarily utilized 
by autonomous kinematically-controlled and driftless^ vehicles [43]. 

^ A kinematically-controlled vehicle is one whose control input is its velocity. A driftless vehicle is a ve¬ 
hicle that "can be stopped instantaneously by setting the control input to zero" [42] - that is, “momentum¬ 
less systems”. Examples of kinematically-controlled and driftless vehicles include autonomous vehi¬ 
cles [38], [43], [44] as well as robotic manipulators under kinematic control [45H48]. 
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Unlike the previous RT-G&C algorithms, Inverse Dynamics in the Virtual Domain (IDVD) 
and Model Predictive Control (MPC) can be considered to be Computational Guidance 
& Control (CG&C) methods [49] as they attempt to solve the G&C task through an iter¬ 
ative, model-based approach. The dividing difference between RT-G&C and CG&C is 
that the latter has the implication of the method being iteration-based, or has otherwise 
been termed as algorithmic control [50] With this being the defining difference, CG&C 
algorithms can be considered a subset of RT-G&C algorithms. The IDVD has been ex¬ 
perimentally evaluated for several tasks, including autonomous rendezvous and dock¬ 
ing (AR&D) and docking with a rotating target [51H53]. The IDVD method is a quasi- 
optimal guidance method that is predicated on parameterizing both the trajectory and 
time as basis polynomials with a set number of coefficients [54]. The resulting poly¬ 
nomial coefficients are then obtained via a nonlinear programming optimization routine 
which minimizes the control effort while satisfying system and path constraints. Due 
to these approximations, the IDVD optimization method can be utilized as a feedback 
controller [52]. In contrast, the MPC framework solves a similar optimization problem, 
but over a limited horizon [55]. This receding horizon control approach allows for the 
efficient computation of the optimal solution to the constrained OCP. Additionally, the 
MPC framework allows for various constraint handling approaches, which can result 
in either a quadratic programming or nonlinear programming problem [52], [56]. The 
MPC approach has been widely applied to real-time spacecraft RPO and rendezvous 
and docking with both a rotating and non-rotating target [56]-[63]. 

Sampling-based methods such as Rapidly-Exploring Random Tree (RRT)-based meth¬ 
ods (described in more detail in Section 6.4.1) offer a hybrid solution which leverages 
the computational simplicity of analytical methods, such APF and AAPF, while achiev¬ 
ing near-optimal performance similar to optimization-based methods, such as IDVD and 
MPC [64], [65]. These methods have been widely used in autonomous kinematically- 
controlled and driftless vehicles, mainly operating in a static environment [44], [64], 
[66]-[68]. Due to their construction, sampling-based methods have the inherent advan¬ 
tage of efficiently exploring the configuration and arriving at a feasible solution relatively 
quickly [64], [67], [68]. 

Lastly, given the vast number of RT-G&C and emerging number of CG&C algorithms 
available in literature, the challenge arises of determining which method is “best” for 
a given application. Specifically, which algorithm enables the system to complete its 
objectives in an efficient manner while simultaneously meeting system-specific require¬ 
ments (e.g., collision avoidance, computational deadlines). Additionally, when reporting 
on a new RT-G&C method, the majority of measurements presented in literature (e.g., 
control effort, maneuvering time, computational time) are implementation dependent. 
These system-specific measurements make it difficult to compare algorithms found in 
literature as the underlying system is fundamentally different (i.e., different physical and 
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actuator properties). Thus, it is apparent that the results be combined to form a metric 
that is agnostic to any specific implementation, denoted here as a “global" metric. The 
resulting global metric can then be used to aid in the selection of candidate algorithms 
to use for a particular application. To the best knowledge of this author, no such global 
metric that characterizes RT-G&C methods has been proposed in literature 


1.3 Research Objectives 

Motivated by the prevalence of RT-G&C algorithms available in literature and the future 
needs of space missions, the main focus of this dissertation is to further the experimen¬ 
tal evaluation methodology for spacecraft guidance and control architectures. This is 
achieved through the following research objectives: 

1. Advancement of a world-class, ground-based, system-level test bed with guar¬ 
anteed real-time capabilities, rigorous characterization, and intuitive software ar¬ 
chitecture for the experimental evaluation and validation of guidance and control 
architectures 

2. Development of a general purpose evaluation framework to compare guidance 
and control methods 

3. Experimental evaluation of the proposed guidance and control method with a 
comparison to those available in literature. 

4. Development of a hybrid guidance and control method based on the Rapidly- 
Exploring Random Tree Star (RRT*) framework that leverages the computational 
simplicity of analytical methods with the optimality of optimization-based meth¬ 
ods. 

The research objectives are cross-referenced with the their supporting chapters in Ta¬ 
ble 1.2. 


1.4 Dissertation Contributions 

The work presented in this dissertation focuses on the advancement of the state of 
the art in the experimental evaluation of Real-Time Guidance & Control (RT-G&C) and 
Computational Guidance & Control (CG&C) methods. This work is fundamental to the 
emerging field of CG&C algorithms as it provides the necessary methodology and tools 
to evaluate their performance. The contributions are summarized as follows: 
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Table 1.2. Mapping of the research objectives to the supporting chapters. 


Research Objective 

Dissertation Chapter 


3 4 

5 6 

7 

Advancement of a world-class, ground-based, 
system-level test bed with guaranteed 
real-time capabilities, rigorous characterization, 
and intuitive software architecture 
for the experimental evaluation and validation 
of guidance and control architectures 

/ / 

/ 

-- 

Development of an general purpose evaluation 
framework to compare guidance and control methods 

-- 

/ / 

/ 

Experimental evaluation of the proposed 
guidance and control method with a 
comparison to those available in literature. 

-- 

/ 

/ 

Development of a hybrid guidance and 
control method based on the RRT* 
framework that leverages the computational 
simplicity of analytical methods with the optimality 
of optimization-based methods. 

/ 

/ 

-- 


• Extensive characterization of the NPS POSEIDYN testbed, a world-class, ground- 
based, system-level GNC test bed. 

• Development and implementation of the TRIDENT-GNC software architecture, 
enabling research into the emerging Computational Guidance & Control field. 

• Development of a standardized evaluation framework to evaluate autonomous 
guidance and control (G&C) algorithms in an equitable manner. 

• Development and validation of a global Guidance Comparison metric that simul¬ 
taneously captures the algorithm performance and computational complexity of 
G&C algorithms. 

• Experimental evaluation of the Sigma-Delta modulator as a thruster modulation 
technique. 

• Application of a proposed RRT-based algorithm as a feedback controller for the 
translational G&C of a spacecraft in a dynamic environment. 

• Development and validation of a metric to estimate the minimum sampling rate of 
a control system. 
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Development and formulation of a G&C trade space for spacecraft design and 
GNC selection. 


The contributions of this dissertation are mapped to the research objectives in Ta¬ 
ble 1.3. 


Table 1.3. Summary of contributions mapped to research objectives. 


Contribution 


Research Objective 


1 

2 

3 

4 

Extensive characterization of the NPS 

POSEIDYN test bed 

/ 

-- 

/ 

-- 

Development and implementation of the 
TRIDENT-GNC software architecture, 
enabling research into the emerging 

Computational Guidance & Control field 

/ 

-- 

/ 

-- 



Development of a standardized evaluation 
framework to evaluate autonomous G&C 
algorithms in an equitable manner 

-- 

/ 

/ 

-- 

Development and validation of a global 

Guidance Comparison metric which 
captures the performance and computational 
complexity of algorithms 

-- 

/ 

/ 

-- 

Experimental evaluation of the Sigma-Delta 
modulator as a thruster modulation technique 

-- 

-- 

-- 

/ 

Application of a RRT-based algorithm 
as a feedback controller for the 

G&C of a spacecraft in a dynamic environment 

-- 

-- 

/ 

/ 

Development and validation of a metric to 
estimate the minimum sampling rate 

-- 

/ 

-- 

/ 

Development and formulation of a G&C 
trade space for spacecraft design and 

GNC selection 

-- 

/ 

/ 

-- 


1.5 Dissertation Overview 

This dissertation is constructed in such an order as to produce a logical flow. After the 
dynamics, constraints, and performance measures associated with a relevant scenario 
are presented, the dissertation details the specifics of the NPS POSEIDYN test bed and 
work performed to increase its robustness, reliability, and repeatability. After character¬ 
izing the test bed and the simulated spacecraft (i.e., test vehicle), the challenge associ¬ 
ated of efficiently modulating the commanded control signal into a sequence of binary 
pulses to control the simulated spacecraft’s thruster is explored. Next, a framework 
is presented which provides an equitable methodology to compare various RT-G&C 
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methods. Several guidance and control methods, including the proposed RRT-based 
method, are then detailed, evaluated, and compared using the proposed evaluation 
framework. The dissertation is then concludes in the traditional manner with a sum¬ 
mary of contributions and recommendations for future work. 

A detailed overview of the remaining chapters is as follows: 

• Chapter 2 introduces the reference scenario, the associated dynamics, relevant 
system and path constraints, as well as the key (optimal) performance measure¬ 
ments. These dynamics, path constraints, and performance measurements will 
be utilized in the development and testing of RT-G&C algorithms. 

• Chapter 3 presents an overview and characterization of the test bed which will 
be utilized to conduct the various experimental campaigns throughout this dis¬ 
sertation. The motivation for ground-based test beds is presented, followed by a 
survey of comparable test beds. A detailed description of the hardware and soft¬ 
ware architectures of the test bed is then presented. Next, the operating system 
latencies, sensor noise, thruster performance, physical properties, and end-to- 
end residual accelerations are characterized. Lastly, a case study is performed 
in order to demonstrate the capabilities of the test bed. 

• Chapter 4 investigates and experimentally evaluates and compares a new thruster 
modulation technique, the Sigma-Delta Modulator (which was briefly considered 
by Ciarcia et al. [69]), to the Pulse-Width Modulator. Linearized models of each 
modulator are presented. Next, these models, in combination with a simple PD 
control law, are used to develop candidate close-loop models. An experimen¬ 
tal campaign is performed using the POSEIDYN test bed, comparing the two 
modulators in a relevant environment with hardware-derived noise and uncer¬ 
tainties. Finally, the results from the experiments, along with various numerical 
simulations, are analyzed to draw conclusions regarding the applicability of the 
candidate models as well as the Sigma-Delta Modulator’s merit as a thruster 
modulation technique. 

• Chapter 5 establishes a proposed experimental framework in which to test and 
(partially) validate candidate guidance and control methods. This framework con¬ 
sists of several relevant test scenarios which attempt to extract various perfor¬ 
mance characteristics and qualities of the candidate algorithm. Next, a unique 
global comparison metric is proposed that captures key performance character¬ 
istics and allows them to be compared simultaneously. 
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• Chapter 6 compares the proposed RRT-based feedback guidance and control 
technique against several other RT-G&C found in literature ranging from analyti¬ 
cal to optimization-based methods. First, each algorithm being compared is intro¬ 
duced and summarized. Next, a detailed overview of sampling-based methods is 
presented, followed by a detailed development of the proposed guidance method. 
The individual results from the experimental campaign are presented. Finally, the 
proposed CG&C method is compared against the other RT-G&C methods evalu¬ 
ated. 

• Chapter 7 provides a case study, using the results from Chapter 6, illustrating 
the use of the Guidance Comparison metric in selecting a G&C method for a 
notional mission. To setup the guidance and control trade space, the Guidance 
Comparison metric is first bounded (in terms of control effort) using parameters 
typically available early in the design process. Next, time-related bounds are 
determined. In particular, a metric predicated on first-principles is introduced to 
estimate the minimum sampling rate required for a given guidance or control task. 

• Chapter 8 concludes the dissertation with a summary of the work performed, 
academic contributions, as well as areas for future work. 


1.6 Publication History 

A complete listing of all conference and journal publications, patents, as well as relevant 
supplemental material generated are listed in Appendix C. Table 1.4 cross-references 
the associated chapter or section to the relevant publication from which content was 
used. 
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Table 1.4. Summary of relevant publications utilized throughout the dis¬ 
sertation. 


Reference 

Number 

Chapter/Section 

Journal/Conference 

Status 

[31] 

Chapter 3 

AIAA JSR^ 

Published 

[33] 

Chapter 3 

2016 AIAA/AAS ASC^ 

Published 

[70] 

Chapter 4 

AIAA JGCD^ 

Accepted 


Section 5.2 



[71] 

Section 6.3.1 

IEEE TCST^ 

Under Review 


Section 6.7 



[32] 

Section 6.3.1 

AAS/AIAA^ 

Published 

[52] 

Section 6.7 

(2016) 6th ICATT® 

Published 

[72] 

Section 6.7 

CEAS Space Journal 

Published 

[61] 

Section 6.7 

AIAA/AAS ASC^ 

Published 

[73] 

Section 6.4 

Section 6.4.2 

2017 AAS/AIAA SFM^ 

To Be Published 

[74] 

Section 7.3.1 
Appendix A 

2017 AAS/AIAA SFM^ 

To Be Published 

[75] 

Chapter 5 
Chapter 7 

AIAA JGCD^ 

To Be Submitted 


^ Journal of Spacecraft and Rockets 
^ Astrodynamics Specialist Conference 
^ Journal of Guidance, Control, and Dynamics 
^ Transactions on Control Systems Technology 
^ Spaceflight Mechanics Conference 

® International Conference on Astrodynamics Tools & Techniques 
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CHAPTER 2: 

Constrained Dynamic System Development 


A problem well put is a problem half solved. 

—Anonymous 
as quoted by Donald E. Kirk 
Optimal Control Theory [ 76 ] 


2.1 Reference Scenario 

The necessity to perform autonomous rendezvous and docking (AR&D) to both a coop¬ 
erative and uncooperative resident space object (RSO) is a fundamental requirement 
for many future missions [3], [4], [46]. Currently, AR&D is routinely performed between 
the ISS using vehicles such as the European Space Agency’s ATV and the Russian 
Soyuz and Progress vehicles. However, vehicles such as the Japan Aerospace Explo¬ 
ration Agency’s H-ll Transfer Vehicle (HTV) and SpaceX’s Dragon capsule are berthed 
to the ISS after autonomously rendezvousing. The reference scenario considered is a 
terminal rendezvous and docking with the ISS with the parameters listed in Table 2.1. 
The remainder of the chapter will detail the dynamical environment and the associated 
equations of motion, relevant system and path constraints associated with AR&D, as 
well as various (optimal) performance measures. 

Table 2.1. Relevant orbital parameters for the reference rendezvous and 

docking scenario Adapted from [77]. 


Parameter 

Value 

Orbital Altitude 

400 km 

Eccentricity 

0.00169 

Inclination 

51.950 deg 

Mean Motion 

0.00113 rad /s 


2.2 Dynamical Environment 

2.2.1 Translational Dynamics 

The relative translational dynamics between two spacecraft, referred to as the Target 
and Chaser (or Chief and Deputy in other literature) spacecraft (denoted in Figure 2.1), 
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is described by the Hill-Clohessy-Wiltshire (HCW) framework. In this framework, it 
is assumed the Target spacecraft has a near-circular orbit. The dextral orthonormal 
basis for the HCW framework is referred to as the RSW frame or Local-Vertical Local- 
Horizontal (LVLH) frame in literature whose axes are anchored at the Target’s center of 
mass (COM) and are defined as follows [78]: 


X ;= 


Z ;= 


Y ;= 


I’target 

Iktargetil2 

I’target ^ V^gi-get 
Iktarg et ^ Vtarget|l2 

Z X X 


Z X X 


(2.1a) 


(2.1b) 


(2.1c) 


where rtarget and vtarget is the respective position and velocity of the Target spacecraft 
in the Earth-Centered Inertial (ECl) frame, defined by the unit vectors (i, J,K)in Fig¬ 
ure 2.1. Note, the unit vectors (i, J,K) of the Earth-Centered Inertial (ECl) frame are 
anchored at the center of the Earth and are defined as follows: the i axis is in the 
direction of the vernal equinox (T); the K axis extends from the origin toward the north 
pole; the J axis completes the dextral orthonormal basis [79]. 


The relative translational dynamics between the Chaser and Target using the HCW 
framework becomes a function of the Target spacecraft’s mean motion, n, and are 
given as [78], 


X = fx + 2ny + 2n^x 


(2.2a) 


y = fy- 2ni; 

z = fz- n^z 



(2.2b) 

(2.2c) 

(2.2d) 


where fx,fy,fz are forces acting on the spacecraft in the x, y, and x directions re¬ 
spectively; /i0 is the gravitational parameter of Earth and a is the semi-major axis of 
the Target spacecraft’s orbit. 


During the terminal rendezvous and docking phase, the Chaser is in close proximity 
of the Target with small relative velocities and the resulting motion is assumed to be 
planar. Additionally, over short durations of time relative to the orbital period of the 
Target, the resulting dynamics can be sufficiently approximated as double integrator 
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T 



Figure 2.1. Illustration of the RSW frame for a Target spacecraft in a near¬ 
circular orbit with the Chaser in close-proximity. Adapted from [80]. 


dynamics [81]: 


x = f. (2.3a) 

y = fy (2.3b) 


2.2.2 Rotational Dynamics 

The attitude dynamics of the Chaser follow Euler’s Equation [82], 

T = Jda (2.4) 

where r G is the net torque acting on the spacecraft; J G is the inertia matrix 
with respect to the COM in the spacecraft body Cartesian coordinate system; lu g 
is the angular velocity of the spacecraft body frame with respect to the inertial frame 
expressed in the spacecraft body Cartesian coordinate system; and is the metrical 
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representation of the cross-product operator, defined as, 


0 

—03 

CI2 

as 

0 

— Cl\ 

— 0-2 

di 

0 


(2.5) 


When the inertia is measured in the Principal Axis Cartesian coordinate system, the 
inertia matrix becomes diagonal. For small angular velocities, the gyroscopic term 
becomes negligible and Equation (2.4) becomes 

T = Jlu (2.6) 

This implies the motion about each axis is decoupled and can be controlled indepen¬ 
dently. Furthermore, assuming no external disturbances, r becomes the control torque. 

For the planar translational motion assumed in Equation (2.3), it is only necessary to 
consider the rotational dynamics about the 2 ; axis. Therefore, the resulting rotational 
dynamics can be simplified to, 

Tz = Jzz^z (2.7) 

where J^z is the inertial about the 2 :-axis in the spacecraft body frame and is the 
angular velocity about the 2 ;-axis of the spacecraft body frame with respect to the inertial 
frame. 


2.3 Constraints 

The constraints on the Chaser spacecraft can be broken up into two main categories: 
dynamical system constraints and path constraints. Constraints on a dynamical system 
(i.e., system constraints) are limitations imposed on the evolution of the system. These 
limitations are typically manifested through actuator saturation limits, such as the min¬ 
imum impulse and maximum force a thruster can generate or the maximum torque 
generated by a momentum exchange device [83]. Path constraints, on the other hand, 
are typically used to describe regions of the state space that are considered admissible. 

2.3.1 System Constraints 

The system constraints considered herein are on the admissible controls. The set 
of all admissible controls, U, considered herein is defined as the set of all controls 
in the inertial frame, ^u(t) g (where Nu is the number of control inputs) such 
that the control input in the spacecraft body-fixed frame, ^u(t), is bounded by some 
maximally negative control input u and some maximally positive control input u+. 
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This admissible set can be written as 


( 2 . 8 ) 


U := {^u{t) e I u- < ®u(t) < u+} 


2.3.2 Path Constraints 

The path constraints considered can be further broken up into two sub-categories: a 
keep-out zone which encompasses an obstacle and a docking cone corridor. These 
two types of path constraints are illustrated in one possible configuration in Figure 2.2. 



Figure 2.2. Notional illustration of a keep-out zone combined with the 
docking cone corridor constraint for the rendezvous and docking scenario. 
Adapted from [84]. 

The first constraint, the obstacle keep-out zone, encompasses the geometry of the 
object and defines an inadmissible portion of the configuration space. In order to satisfy 
the constraint, the Chaser COM must remain outside of the obstacle keep-out zone. 
The keep-out zone constraint can be formally stated as 

ll^ci||2 — (^ci) (2-9) 

where Vd is the relative position of the Chaser COM with respect to the COM if the 
obstacle and di is the distance from the COM of the obstacle to the boundary of its 
keep-out zone along the line r^. Note, as indicated by Equation (2.9), di is not constant 
and allowed to vary as a function of the relative position of the Chaser with respect to 
the obstacle. The only requirement placed on di is that the surface it describes is 
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continuous. One such instantiation of a keep-out zone is illustrated Figure 2.2. Other 
examples of keep-out zones are illustrated in [19], [32], [85], [86]. 

The second type of path constraint considered is a docking cone corridor (DCC). This 
path constraint ensures a safe, collision-free final trajectory to the Target in addition to 
ensuring the docking action occurs parallel to the docking axis [28], [87]. The DCC 
is anchored at the docking interface of the Target spacecraft and is defined by three 
parameters: the docking axis unit vector, riDcc, the docking cone half-angle, (3, and the 
length of the corridor, 4- The DCC path constraint is given as 

||rcd||2COs(/3) < nJccrcd (2.10) 

where r^d is the relative position of the Chaser with respect to the Target docking in¬ 
terface. Unlike other constraints, the DCC constraint is only active when the distance 
between the Chaser COM and the Target docking interface is less than or equal to 4 - 


2.4 Optimal Control Problem Formulation 

Optimal control is predicated on determining the control input u g U C necessary 
to transfer a system x g from an initial state to a desired final state that minimizes a 
specified cost functional [88]. Through the application of the Calculus of Variations, the 
necessary conditions can be derived through the introduction of the costates [76], [88]. 
However, the dimensionality of the resulting optimization problem is increased by a fac¬ 
tor of two from to 2Nx. Notably, this increase in dimensionality is commonly referred 
to as the Curse of Dimensionality [76], [88]. While the increased dimensionality is not 
necessarily an issue for modern (terrestrial) computational technology, implementation 
onboard a spacecraft provides unique challenges due to the computational constraints 
of space-based compute elements. In the remainder of this section, the general op¬ 
timal control problem will be developed, along with two common optimal performance 
measures: the minimum-propellant and minimum-time optimal control problems. 

The general optimal control problem can be succinctly stated as the following [76], [88], 
[89]: 


r*f 

Minimize: J(-) = /i (x(t/),t/) + / g {x{t), u{t), t) dt 

Jto 

(2.11a) 

Subject To: X = f (x(t), u(t), t) 

(2.11b) 

x(to) = X° 

(2.11c) 

Cp (x(t),u(t),t) < 0 

(2.1 Id) 

e/(x/,t/) = 0 

(2.11e) 


18 



where the cost functional J(-) : x h-)- M is composed of the endpoint cost 

h{-) : R^^ X M !-)■ M and the running cost g{-) : R^^ x x M i-)- M. The system 
dynamics x of the system is given by f : R^^ x R^^R i-)- R^^ with an initial condition 
x° at time to and terminal condition ej : R^^ x M i-)- Lastly, any constraints, such 
as path and/or actuator constraints, are given as Cp : R^^ x x M i-)- 

Using the Calculus of Variations, the necessary conditions for optimality for a system 
with bounded control inputs, u g U, can be derived. To begin, the cost functional given 
in Equation (2.1 la) can be rewritten as an augmented cost funcitonal Ja(-) to include 
differential constraints [76], 


■U-) = 


rtf 


'to 



dh , , 

!//(■) + 



f)h 

Mt) + -^i-) + A^(f(-) - x(t)) ) dt 


( 2 . 12 ) 


rtf 


9a{-)dt 


'to 


where \{t) e R^^ are the costates (which are also referred to as Lagrange Multipliers 
in some literature). Defining the Hamiltonian as [76], [88], 


,^(x, u, A, t) = g{-) + A(y/(■) 


(2.13) 


the augmented integrand ga{-) becomes 


gai-) = ^i-) - A^x(t) + 


dh 


) x(/) + 



(2.14) 


Invoking the Fundamental Theorem of the Calculus of Variations [76], the necessary 
condition under which the augmented cost functional, with the augmented integrand 
given by Equation (2.14), is minimized is given as [76] 

6Ja{u*,6u)>0 V ugU (2.15) 

where u* is the optimal controP that minimizes the cost functional J(-). Note, for the 
case of an unbounded control, the inequality relation in Equation (2.15) becomes an 
equality relation 

(5Ja(u*,5u) = 0 V UGU 

''Components associated with the optimal solution are denoted with an asterisk in the superscript (•)* 
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since the u* lies within the control set bounded by ±cx) [76]. Evaluating the variation of 
the cost functional in Equation (2.15), the necessary conditions for optimality are [76], 


= aA <■> 

(2.16a) 

ax<'> 

(2.16b) 

u*, u, 

(2.16c) 

with the boundary condition given as 


dh 



dh ^ 

+ 5x^ = 0 (2.17) 


where 5tf is the variation of the final time and (5x/ is the variation of the terminal state. It 
is important to note. Equation (2.16b) specifies the costate dynamics. Equation (2.16c) 
is known as Pontryagin’s Minimum Principle, and when the boundary conditions are 
specified. Equation (2.17) reduces to the (terminal) transversality condition [76], [88]. 

2.4.1 Minimum-Propellant Optimal Control Problem 

The minimum-propellant optimal control problem is arguably one of the most prevalent 
OCPs for RPO. As its moniker implies, the resulting solution from solving this OCP is 
the one that uses the least amount of fuel. Similar to the development of the generic 
optimal control problem, it is first assumed that the control input u g U C is 
bounded. The cost functional associated with the minimum-propellant optimal control 
problem is the Li-norm of the control input, ||u(t)||^ and can be defined as [76], 

Nu 

l|u(t)ili = 5^|Mi(t)| (2-18) 

i=l 

The generic optimal control problem given in Equation (2.11) can be rewritten to form 
the minimum-propellant optimal control problem: 

r^f 

Minimize: J(-) = / ||u(t)||^ dt (2.19a) 

Jto 
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Subject To: X = f (x(t), u(t), t) (2.19b) 

x(to) = x° (2.19c) 

Cp (x(t), u(t), t) < 0 (2.19d) 

e/(x/,t/) = 0 (2.19e) 

For the purpose of illustration, let the system dynamics be defined as [76], 

X = a(x(t),t) + B(x(t),t)u(t) (2.20) 


where 

B(-) = lbi(-),b2(-),...,bjv^(-)] 

may depend on the states and time.^ Note, the assumption of the system dynamics 
being affine in control is common for dynamical systems [76], [90]. The resulting Hamil¬ 
tonian (using Equation (2.13)) associated with the minimum-propellant optimal control 
problem is 

e^(x, u, A, t) = ||u(t)||^-|-A(t)'''a(x(t), t)-|-A(t)'''b(x(t), t)u(t) (2.21) 


Applying Pontryagin’s Minimum Principle (given by Equation (2.16c)) to the minimum- 
propellant OCP yields 

||u*(t)||^ -h A*(t)''’a(x*(t),t) -h A*(t)’'’B(x*(t),t)u*(t) < 

l|u(t)||i + A*(t)Vx*(t),t) -F A*(t)^B(x*(t),t)u(t) (2.22) 

Assuming B(-) is diagonal (implying no interactions between the components of u) 
Equation (2.22) can be simplified to [76], 

K {t) I + u* {t)X* (t)^bi(x* (t) , t) 

<\u,{t)\+Ui{t)X*{tybi{x*{t),t) , t = l,2,...,Nu (2.23) 

‘'When B(-) is not dependent upon any state or time, it reduces to the control input matrix traditionally 
associated with linear time-invariant systems. 
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The conditions under which the Hamiltonian is minimized (and Equation (2.23) is satis¬ 
fied) can be ascertained by expanding the right-hand side of Equation (2.23) [76]: 


\ui(t)\ + Ui(t)X*(tyhi(x*(t),t) 


[l + A*(t)'''bi(x*(t), t)] Mj(t) forMj(t) >0 
[-1 + X*{tyhi{x*{t),t)]ui{t) ior Ui{t) <0 


(2.24) 


By inspection, the optimal control law can be found to be [76]: 



ut 

for 

x*{tyh,{^*{t),t) < -1 

Ui 

for 

A*(t)’'’bi(x*(t),t) > 1 

0 

for 

x*{tyh,{^*{t),t) E (-1,1) 

undetermined 

for 

A*(t)Tb,(x*(t),t) = ±1 


(2.25) 


where M+ is a non-negative number (i.e., M e [0, cx))) and M<o is a strictly negative 
number (i.e., M e (-cx),0)). It is worthwhile to note, an interval is said to be singular 
when X*(tyh{-x*(t),t) = ±1 for a nonzero duration of time; whereas a crossing of ±1 
merely indicates a control switch from u~ —)■ or vice versa [76]. The resulting control 
action specified by the optimal control policy in Equation (2.25) is termed bang-bang or 
bang-off-bang control, assuming no singular intervals. 

2.4.2 Minimum-Time Optimal Control Problem 

Given a bounded control input u g U C the minimum-time optimal control problem 
is given as: 

r^f 

Minimize: J(-) = / Idt (2.26a) 

Jto 

Subject To: X = f (x(t), u(t), t) (2.26b) 

x(to) = x° (2.26c) 

Cp (x(t), u(t), t) < 0 (2.26d) 

= 0 (2.26e) 

Again, for the purpose of illustration, let the system dynamics x be 

X = a(x(t), t) + B(x(t), t)u(t) (2.27) 

Note, unlike for the minimum-propellant OCR, the minimum-time OCR requires the final 
time be free. From Equation (2.13), the Hamiltonian is associated with the minimum- 
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time OCP is 


e^(x, u, A, t) = 1 + A(t)'''a(x(t), t) + A(t)'''b(x(t), t)u(t) (2.28) 


Applying Pontryagin’s Minimum Principle to the minimum-time OCP yields 

1 + A*(t)'''a(x*(t),t) + A*(t)'''B(x*(t),t)u*(t) < 

1 + A*(t)'''a(x*(t), t) + A*(t)'''B(x*(t), t)u(t) 


which can be simplified to 

i = l,2,...,Nu (2.29) 

Following a similar procedure to determine the minimum-propellant optimal control pol¬ 
icy, the optimal control policy for the minimum-time problem is 

( ul for A*(t)’'’bi(x*(t),t) < -0 

u*{t) = < u- for X*ltyhilx*lt),t) > 0 (2.30) 

[ undetermined for A*(t)^bj(x*(t), t) = 0 

It is worthwhile to note, an interval is said to be singular when A*(t)^b(x*(t), t) = 0 for 
a nonzero duration of time; whereas a crossing of ±1 merely indicates a control switch 
from u~ -)■ or vice versa [76]. 

It is worthwhile to highlight another property of time-free OCPs by examining the bound¬ 
ary conditions given in Equation (2.17). Since the terminal state is specified, the bound¬ 
ary conditions for minimum-time OCP is given as 

= o (2.31) 

Since Stf y 0 since the final time is not specified, = 0 in order to satisfy the 

boundary condition. Furthermore, since the Hamiltonian is not an explicit function of 
time, its time-derivative is zero which implies 

,^(x*(t),u*(t),A*(t),t) = 0 V te[to,tf] (2.32) 

This fact that the Hamiltonian is zero when evaluated on the extremal trajectory can be 
used in validating the optimal solution [76]. 
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CHAPTER 3: 

The POSEIDYN Test Bed and TRIDENT-GNC 
Software Architecture 


I constructed a laboratory in the neighborhood of Pike’s Peak. The 
conditions in the pure air of the Colorado Mountains proved 
extremely favorable for my experiments, and the results were 
most gratifying to me. 

—Nikola Tesla 
“Talking with the planets” [91] 


The work presented in this Chapter was originally presented at the American Institute 
of Aeronautics and Astronautics (AIAA)/American Astronautical Society (AAS) Astrody- 
namics Specialist Conference 2016 in Long Beach, CA (12-16 September 2016) [33]; 
a journal version of this work has been published as an article in advance in the AIAA 
Journal of Spacecraft and Rockets [31 ]. 


3.1 Motivation for Ground-Based Test Bed 

Integration, verification, and validation are some of the greatest challenges which fu¬ 
ture autonomous GNC systems face [2], [3]. To ensure safe on-orbit operations for 
RPO, ground testing must be performed to ensure the resulting GNC algorithms meet 
predefined set of performance requirements. To perform this task, at the recommen¬ 
dation of Quadrelli et al., it is argued that ground-based, system-level, demonstration 
systems can provide adequate fidelity to verify and validate GNC performance require¬ 
ments [4], Air-bearing test beds provide a dynamically representative environment to 
develop and test GNC algorithms. In these facilities, the test vehicles that represent 
spacecraft (or resident space objects) operate on top of a planar surface. Air-bearings 
located on the test vehicles are used to reduce the friction of the vehicles to a cre¬ 
ate a quasi-frictionless environment. As the planar surface is horizontally leveled, the 
effects of any in-plane components of gravity on the test vehicles are minimal. The 
result is a quasi-frictionless and low residual acceleration environment in a plane, thus 
emulating the drag-free and weightless environment of orbital spaceflight. Additionally, 
these test beds include various hardware phenomena such as delays, computational 
constraints, actuator response uncertainty, and sensor noise, which is practically im¬ 
possible to replicate in a simulated environment. Ground-based GNC test beds are a 
useful tool to advance the state-of-the-art of these systems and can also be used to 
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perform end-to-end system-level verification and validation prior to the system’s opera¬ 
tional deployment. 


3.2 Survey of Relevant Test Beds 

Although survey papers on air bearing spacecraft simulators for RPO exist in literature 
[48], [92], new facilities have emerged since these were published. A summary of the 
characteristics of the examined test beds not found in the survey literature is tabulated 
in Table 3.1. 


Table 3.1. Summary of the relevant test beds characteristics. Adapted 
from [93]-[97]. 


Name 

Location 

Simulator Type 

Surface Material 

Navigation Sensors 

ADAMUS 

UF 

Translational + Rotational 
Air-Beaing 

Epoxy 

PhaseSpace 

ASTROS 

Georgia Tech 

Translational + Rotational 
Air-Bearing 

Epoxy 

Vicon 

SICK Laser sensor 
IMU 

ORION 

FIT 

Translational + Rotational 
Air-Bearing 

Epoxy 

OptiTrack 

PINOCCHIO 

Sapienza 

Translational Air-Bearing 

Glass 

On-board IMU 

POSEIDYN 

NPS 

Translational Air-Bearing 

Granite 

Vicon 

FOG 


Both the University of Florida (UF) ADAMUS test bed and the Florida Institute of Tech¬ 
nology (FIT) ORION (with its DAWN vehicle) test bed are capable of the full six degree 
of freedom (DOF) motion [93]-[95]. In these test beds, planar air bearings are used 
by the test vehicles to translate on top of an epoxy floor [2-DOF] while a spherical air¬ 
bearing provides the rotational 3-DOF and a complex counter-balance system is used 
to provide frictionless out-of-plane motion. The Georgia Institute of Technology AS¬ 
TROS utilizes a 5-DOF platform, where the ASTROS platform is a 2-DOF rotational air 
bearing atop a 3-DOF platform that employing translation air bearings [96]. Lastly, the 
University of Rome-La Sapienza PINOCCHIO test bed is a two-translation and one- 
rotational DOF test bed utilizing only a translational planar air-bearing [97]. 

With the exception of the PINOCCHIO test bed at University of Rome-La Sapienza, 
which is still in development, all the examined test beds utilize a commercial motion 
capture system, such as Vicon or PhaseSpace, in order to provide an inertial position 
to the spacecraft simulator. These external position measurements are usually aug¬ 
mented by an on-board inertial measurement unit (IMU) and fused together for more 
accurate state estimations. 
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Furthermore, the use of an epoxy surface is pervasive throughout the surveyed air¬ 
bearing test beds. Although granite surfaces provide unmatched surface planarity and 
smoothness, epoxy floors or glass surfaces have a significantly lower procurement 
cost, thus making them more common. Lastly, it is worthwhile to note that a detailed 
characterization of the test beds is typically not found in literature, thus hindering the 
comparison of the experimental results obtained when using different test beds. From 
an experimental stand point, having such a detailed description and characterization 
is important in order to understand the system’s performance (e.g., navigation per¬ 
formance) and uncertainties (e.g., actuator performance, moment of inertia estimate, 
residual accelerations) under which a GNC algorithm has been experimentally evalu¬ 
ated. 


3.3 Experimental Facility Overview: 

The POSEIDYN test bed 

The Naval Postgraduate School (NPS) Proximity Operation of Spacecraft: Experimen¬ 
tal hardware-ln-the-loop DYNamic simulator (POSEIDYN) test bed has been developed 
with the aim to provide a dynamically representative system-level platform upon which 
to develop, experimentally test, and partially validate GNC algorithms for RPO. Before 
the installation of the 4-by-4 meter granite monolith in early 2012, the vehicles oper¬ 
ated over an epoxy floor since the test bed was initially stood-up in 2004. Over the 
decade that the test bed has been in operation and in constant upgrade, four different 
generations of floating spacecraft simulators (FSSs) have been developed. The evolu¬ 
tion of these test vehicles is shown in Figure 3.1. Each FSS generation has included 
a unique capability. The first generation FSS utilized a prototype capture system in 
order to perform rendezvous and docking [22]. The second generation FSSs featured 
vectorable thrusters and a miniature control moment gyroscope (CMG) [98], [99]. The 
third generation FSSs moved away from an aluminum construction in favor of a more 
lightweight, polycarbonate structure (fabricated using additive manufacturing) and com¬ 
ponents, such as docking interfaces. Lastly, the fourth-generation FSSs continued the 
use of polycarbonate components and includes a standardized interface for use with 
robotic manipulator research [47]. 


3.4 Hardware Architecture Overview 

The POSEIDYN test bed is composed of three main elements: a 4 x 4 m granite mono¬ 
lith, multiple FSSs, and a laboratory-wide metrology system. An overview of these 
elements is shown in Fig. 3.2. 
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Figure 3.1. Lineup of the first- to fourth-generation FSSs used on the 
POSEIDYN test bed. 



Figure 3.2. Overview of the main elements of the POSEIDYN test bed. 


The FSSs are custom-designed vehicles that emulate orbital spacecraft moving in close 
proximity of another vehicle or object (e.g., another FSS). In this paper, the fourth- 
generation FSSs will be considered. Three flat 25 mm diameter air-bearings are used 
by the floating spacecraft simulator to achieve quasi-frictionless motion on top of the 
granite monolith. The air-bearings use compressed air to lift the FSSs approximately 
5/im, creating an air film between the vehicle and the granite surface that acts as a 
lubricant layer and eliminates their direct contact. The compressed air, supplied at 
4.1 bars (60 psi), is delivered from an onboard tank holding 1.87 L of compressed air 
at 206.8 bars (3000 psi). With a nominal air consumption of 0.1 L/min at 4.1 bars 
(0.53 normalized liters per minute) per air-bearing, the floating endurance of the FSS is 
estimated to be at approximately 3.5 h. A solenoid valve controls the flow toward the air- 
bearings and an air filter prevents contaminants and foreign materials from damaging 
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the delicate air bearing porous material. A scuba-diving-type compressor is used to 
refill the onboard tank. 

The 15 metric ton, 4-by-4 meter granite monolith is supported by three adjustable 
pedestals and has a planar accuracy of ±0.0127 mm. This allows the granite sur¬ 
face to be horizontally leveled to an accuracy of 0.01 deg. Seismic activity and other 
building structural activity can distort the monolith alignment over time. Periodic checks 
are performed to ensure that the granite monolith is leveled within the pre-specified 
tolerance. 

The combination of the reduced friction (provided by the air-bearings) with the low resid¬ 
ual acceleration (provided by the horizontality and planarity of the granite monolith) 
allows the POSEIDYN test bed to recreate the drag-free and weightless environment 
experienced by spacecraft in close proximity maneuvering. A characterization of the 
residual acceleration experienced by the FSS can be found in Section 3.6. An impor¬ 
tant limitation of the test bed consists in allowing only 3-DOF motion (two translation 
and one rotation) compared to 6-DOF motion of a spacecraft on orbit. Furthermore, the 
applicability of the test bed, as a high-fidelity dynamic simulator, is limited to short du¬ 
ration close proximity operations because the relative orbital mechanics are not recre¬ 
ated. Despite these limitations, the POSEIDYN test bed is critical to experimentally 
test proximity operations guidance, navigation and control methods in a realistic envi¬ 
ronment, including sensor and actuator dynamics, as well as real-time computational 
constraints. 

To propel each FSS over the granite monolith the vehicles are equipped with eight cold- 
gas thrusters, with two thrusters mounted equidistant from the center. Each thruster is 
composed of a solenoid valve and a custom-made convergent-divergent nozzle [100]. 
Fed by the on-board tank of compressed air, the thrusters provide a nominal thrust 
between 0.1-0.15 N of thrust (see Section 3.6 for thruster characterization details). 
This nominal thrust fluctuates considerably, as the thrust is a function of nozzle inlet 
pressure, which changes depending on the number of thrusters that are being fired 
simultaneously. With a nominal mass flow of about 0.3 g/s per thruster, the total firing 
time is estimated to be 20 min. Red light-emitting diodes paired with each thruster 
provide a visual indication of which thruster is firing. 

Furthermore, one of the fourth generation FSSs has a Ball Aerospace 2.5 N-m s reaction 
wheel (RW) (model RW 2.5-A1) mounted atop its structure, as illustrated in Figure 3.3. 
This particular FSS is used for spacecraft robotics research, and a multilink modular 
robotic manipulator can be connected to it. The RW was installed to meet the increased 
torque requirements due to the dynamic coupling (arising from the manipulator motion) 
and for the significantly higher inertia of the combined system [47]. An Arduino Due 
microcontroller is used as an interface between the on-board computer and the RW. 
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This microcontroller converts the output of the three hall sensors on the RW to a ro¬ 
tation rate and provides a suitable analog signal to control the RW motor torque. A 
switching shunt regulator is used to safely dissipate the back-EMF generated during 
RW deceleration. 



Figure 3.3. 2.5 N m s Ball Aerospace reaction wheel assembly mounted 
atop one of a FSS. 


The onboard computational capabilities of the FSS are provided by a commercial off- 
the-shelf PC-104 form-factor onboard computer. Based on an Intel Atom 1.6 GFIz 32- 
bit processor, the computer has 2 GB of RAM and an 8 GB solid-state drive. Despite 
the FSS onboard computer not being a space-qualified compute system (such as the 
BAE RAD750 or Proton 200k compute elements), its computational capabilities could 
be regarded to be on the same order of magnitude, in terms of the estimated MIPS- 
per-watt,^ as state-of-the-art space-grade computers, as illustrated in Table 3.2 [101]- 
[107]. An expansion board with 20 optoisolated relays provides the required switching 
capability for the thrusters and air bearing solenoid valves. A serial interface is used 
to communicate with an on-board KVH DSP-3000 fiber-optic gyroscope (FOG), which 
provides angular velocity measurements at a 100 Hz rate. On the FSS with the RW, 
another serial port is used to communicate with the microcontroller interfacing with the 
RW. Two 95 W-h lithium-ion batteries and a battery management module regulate the 
electrical power to the FSS. Under idle conditions, the batteries can power the FSS for 
over 5.5 hours; under maximum loading, including the RW, the FSS can operate for 
just under 3.5 hours, thus making the amount of air inside the onboard tank the limiting 
factor during experiments. 

^ MIPS: Millions of Instructions Per Second 
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Table 3.2. Comparison of space-qualified and typical terrestrial proces¬ 
sors. Adapted from [101]-[107] 



Name 

No. of 

Processor Freq. 

Peak Power 

Est 

MIPS per 


Cores/Threads 

(MHz) 

(W) 

MIPS^ 

Watt 


BAE RAD750 

1/1 

200 

12.5 

400 

32 

Space 

Proton200k DSP 

1/1 

200 

5 

4,000 

800 

Qualified 

Proton200k Lite 

1/1 

66 

1.5 

1,200 

800 

Processors 

Tyvak Intrepid 
(ARM Processor) 

1/1 

400 

0.3 

440 

1,467 


ARM Cortex A15 
(Tegra Processor) 

4/4 

1,900 

2.25 

4,000 

1,778 

Terrestrial Intel Core 17-471OMQ 
Processors (Generic Laptop CPU) 

4/8 

2,500 

47 

40,000 

851 


Intel Atom Processor 
(FSS Processor) 

1/2 

1,600 

5 

3,200 

640 


^ MIPS: Millions ol Instructions Per Second 


VICON 

Server 





Figure 3.4. Overview of POSEIDYN communications architecture. 


A Wi-Fi module provides wireless communication capabilities to each FSS [108]. The 
POSEIDYN test bed general communications arrangement is shown in Figure 3.4. The 
Wi-Fi module enables the FSS to communicate with other FSSs or other external com¬ 
puters, such as the laboratory motion-capture system (Vicon) server and ground station 
for telemetry logging and visualization purposes. To minimize latency the data is trans¬ 
mitted from node to node using User Datagram Protocol (UDP). 

A carbon fiber reinforced polymer base plate and four aluminum T-slotted columns 
make up the core of the FSS structure. A polycarbonate outer shell provides the at- 
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tachment points for all the FSS components. Fabricated using additive manufacturing, 
the outer shell can be quickly modified to accommodate any component changes or 
vehicle upgrades (e.g., this was done with the RW). 

The laboratory motion-capture system is a commercial system (Vicon) composed of 
10 overhead cameras and an external computer. This system determines the position 
of objects carrying passive markers (i.e. the FSS) with sub-millimeter level (static) 
accuracy at rates up to 100 Hz. Once the location of the FSS is determined by the Vicon 
system, an external computer streams the data to the FSS using the Wi-Fi link via UDP. 


3.5 Software Architecture Overview 

At the core of the FSS software architecture is a real-time operating system (RTOS) 
which ensures the overlaying GNC software running on-board responds to sensor in¬ 
puts and generates the appropriate actuator outputs within a strict and predefined time 
span. 

To achieve the desired real-time requirement, a Ubuntu 10.04, 32-bit server edition 
operating system (OS) has been chosen and its Linux kernel (v2.6.33) has been patched 
with the RT-Preempt patch [109], [110]. This particular OS combination used by the 
FSS will be referred as RT-Linux OS. Running atop the RT-Linux OS is the multi-rate 
GNC software, which can be broken-up into four main subsystems: navigation, guid¬ 
ance, control, and telemetry. 

The navigation subsystem provides a full state estimate of the FSS (position, orienta¬ 
tion, and rates). Because most of the research conducted on the POSEIDYN test bed 
is focused on the guidance and control aspect, the navigation subsystem makes use 
of the accurate Vicon information to provide an accurate state estimate. As illustrated 
in Figure 3.5, the navigation block first samples the on-board sensors and the actuator 
states and fuses them via a Discrete Kalman Filter (DKF) to produce an inertial state 
estimate. The inertial state estimate is then sent to the guidance subsystem, where 
the appropriate actuator commands are generated and sent to the control subsystem. 
The standard control subsystem is comprised of the steering logic and converts the 
guidance-desired actuator inputs to the required low level signals to drive the differ¬ 
ent on-board actuators. Typically, only the guidance subsystem is user-defined and 
the navigation and control subsystems utilize predefined software. However, when re¬ 
search is performed on relative navigation, for example, the navigation block can be 
used to provide a "ground truth" estimate to benchmark the experimental relative navi¬ 
gation results [22]. 

Lastly, the telemetry subsystem packages the requested telemetry and sends it to a 
desired ground station. In addition to the GNC telemetry specified in the model, a 
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Figure 3.5. Overview of the onboard software architecture of the Floating 
Spacecraft Simulators. 


publicly available system monitoring utility is used to collect metrics for various system- 
level metrics such as, but not limited to, processor usage, memory consumption, and 
network bandwidth. In particular, the sysstat utility was chosen [111]. It is worthwhile 
to note, the sampling rate of each subsystem (navigation, guidance and control) is 
user-specified, allowing for a multi-rate GNC formulation. 

The remainder of this section will detail specifics about the GNC software development, 
navigation formulation, thruster mapping and modulation, reaction wheel control, and 
the development simulator. 

3.5.1 TRIDENT-GNC Development Environment 

To simplify the algorithm development and subsequent implementation on the FSSs, a 
(numerical) development simulator and a FSS software template were created using a 
common custom library. This library contains common software (navigation and control 
subsystems) which is used in both the simulator and the FSS autogenerated onboard 
software. The simulator uses simulated sensors and actuators and also simulates the 
plant response (i.e., the FSS); whereas the FSS software template uses the interfaces 
to the onboard sensors and actuators. This commonality between the simulator and 
FSS onboard software allows for rapid development of the algorithms in a simulation 
environment and software generation for use onboard the FSS for testing. 

The multi-rate GNC software running atop the RT-Linux OS is developed utilizing the 
MATLAB and Simulink R2015b environment. Once developed, the Simulink models are 
autocoded into C and compiled (using the ert.linux target language compiler [112]). 
To facilitate the code generation and compilation across multiple OSs and architec- 
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tures, an Ubuntu virtual machine (VM) has been created the development tool-chain 
pre-installed. This VM allows the user to easily create hard real-time, multi-rate GNC 
software. The aggregation of the custom library, (numerical) simulator, FSS software 
template. Standard Test Framework test scenarios, and development VM and toolchain 
are collectively referred to as the Toolkit for Real-time DevelopmENT of Guidance, Nav¬ 
igation & Control (TRIDENT-GNC) (whose logo is presented in Figure 3.6. 



POSEIDYN 


Figure 3.6. TRIDENT-GNC logo. 


3.5.2 Navigation 

A DKF fuses the Vicon and FOG data as well as an estimate of the actuated force and 
torque (given the states of the thruster valves and RW torque) to provide an estimate 
of the FSS inertial state. The FSS DKF is conceptually broken up into an outer-loop 
consisting of the prediction steps and an inner-loop composed of the correction steps 
that update the outer-loop while the filter is converged. Furthermore, in order to re¬ 
tain flexibility, the DKF operates at a user-defined time step, AtoKF- However, dropped 
or corrupted sensor measurements are inevitable due to the asynchronous sensors, 
user-defined time step, and lack of error correction in UDP communications. A dy¬ 
namic construction of the measurement mapping matrix is used to mitigate the effect 
of dropped sensors measurement and a x^-Qating rejects corrupted out-of-bound mea¬ 
surements [113]. Lastly, in addition to the standard output of a state estimate and error 
covariance estimate, additional telemetry indicating the status and “health” of the filter 
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(e.g., initialization status, convergence status, measurement validity, and measurement 
gating) is stored in the bit field of a 32-bit unsigned integer, illustrated in Figure 3.7. This 
section will detail the construction and considerations given in the development of the 
FSS DKF. 



(MSB) (LSB) 


MO: No Vicon Measurement 
M1: No FOG Measurement 
GO: X Measurement Gated 
G1: y Measurement Gated 
G2: 9 Measurement Gated 
G3: 9 Measurement Gated 
SO: DKF Initialized 
S1: DKF Covariance Converged 

Figure 3.7. Description of the DKF health telemetry item bit fields. 

To formulate the DKF, the system is assumed to take the discrete-time representation 

Xfc = + Ffc-iUfe-i -|- (3.1) 

r • 1 

where the state vector is denoted by x*. = 9,9k e at sample time 

k\ Ufc_i = [/a;,/j/,r]^ e is the control input at the previous sample time k - l\ 
^k e is the state transition matrix; r^-i G is the discrete-time input 

gain matrix; g is the process noise covariance matrix; and ujk G R^^ 

denotes the process noise which is assumed to be zero-mean Gaussian white noise. 
Additionally, measurements acquired at sample time k, Zk G R^^, are assumed to 
follow 

Zk = HfcXfc + Vk (3.2) 

where G is the measurement mapping matrix and Vk G is the sen¬ 

sor noise that is assumed to be zero-mean Gaussian white noise. For this particular 
application, the overhead motion capture system (Vicon) provides position and attitude 
measurements while the onboard FOG provides attitude velocity information. 
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The predicted state estimate, ^k\k-i, and predicted error covariance, Pfc|fc-i, are given 
as follows: 


^k\k—l 1 “h PfcUfc—1 (2'2) 

Pfc|fc-i = + Qfc (3.4) 


where the process noise covariance e used in the propagation of the 

predicted error covariance is assumed to have the following form [113], [114], 


Qfc 


^^dkf/4 ^^dkf/3 0 0 0 0 

^^dkf/3 ^^dkf 0 0 0 0 

0 0 AtQ(^p/4 Atpi^p/2 0 0 

0 0 AtQ(^p/2 Atpi^p 0 0 

0 0 0 0 AtQi^p/4 AtQ(^p/2 

0 0 0 0 AtQi^p/2 ^^dkf 


(3.5) 


For simplicity, it is assumed the output of the thrusters follow a square wave - that is, 
the output of the thrusters is modeled to be fully ON or fully OFF with no transient re¬ 
sponse. In order to generalize the implementation of the filter, a forward initialization 
method is utilized to provide an in-situ initial state estimate and initial error covariance 
using two complete measurement sets [113]. Furthermore, to compensate for uncer¬ 
tainties associated with thrusting (see Section 3.6), the noise intensity, a strictly positive 
quantity, ql e M+ associated with the attitude rate process noise diagonal element [see 
Equation (3.5) ] is increased by a factor of 500 from a nominal value of 1 x 10“^. This is 
done because the attitude rate, due to the small moment of inertia of the FSS, is more 
sensitive to thrust uncertainties. Note, the inflation factor was chosen through testing 
on the hardware. 


The next consideration in the construction of the FSS DKF is the dynamic development 
of the measurement mapping matrix g Nominally, the measurement 

mapping matrix is defined as follows: 




1 0 0 0 0 0 
0 0 1 0 0 0 
0 0 0 0 1 0 
0 0 0 0 0 1 


(3.6) 


The first three rows of the Equation (3.6) are associated with measurements from the 
Vicon sensor, whereas the last row is associated with measurements from the FOG. 

The dynamic construction of measurement mapping matrix involves ensuring each 
measurement received from each sensor is “valid” before it is processed. The data 
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validity check is a two-step process; first, the Data.Available flag from the Vicon sen¬ 
sor and FOG are checked to be true. If either (or both) Data.Available flags are 
false, the respective rows in the nominal measurement mapping matrix in Equa¬ 
tion (3.6) are set to a row of zeros and the respective bits in the DKF health telemetry 
item are set to 0. It is important to note that the Data.Available flags are an output 
from the implemented sensor blocks. Measurements which pass the first data validity 
check are then subjected to a measurement association check in order to protect 
against corrupted measurements and filter smugness - a phenomenon which occurs 
when the estimated error covariance shrinks causing the filter to effectively reject new 
measurements resulting in a diverging state estimate unless corrected [78]. Therefore, 
a given measurement is considered to be associated with the FSS to a user-defined 
probability when [113], 

(zfc ^ (3.7) 

where C is a distributed random variable with one degree of freedom and g 
]^NmxNm jg t |-|0 sensor noise covariance matrix and is defined as follows: 


Rfc 


al 0 0 O' 

0 0 0 

0 0 0 
0 0 0 

u - 


(3.8) 


Note, since each measurement is assumed to be independent, there is only one degree 
of freedom (DoF). Additionally, the default value chosen is 99.95%, which results in 
C 99 . 95 % = 12.116. Although Equation (3.7) is written in metrical form for compactness, 
it can be easily expanded into four scalar equations for easier comparison with the 
association threshold. Since the measurement matrix does not contain any off-diagonal 
terms, the computations for each measurement are independent, thereby eliminating 
costly matrix inverses. If any measurement violates Equation (3.7), that specific row in 
the measurement mapping matrix is set to all zeros and the respective bit in the DKF 
health telemetry item for gating is set to 1 . Additionally, in an effort to overcome filter 
smugness, the diagonal term in the error covariance estimate is artificially inflated by a 
user-defined scalar - whose default value was chosen through testing to be 1.025 - in 
a similar fashion to a Fading-Memory Kalman Filter [114], [115]. 

Once the measurement mapping matrix is created, the state estimate correction and 
error covariance correction are computed in the traditional manner where the Kalman 
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gain, state estimate, and error covariance estimate are. 


( 3 . 9 ) 

(3.10) 

( 3 . 11 ) 


Kfc — P(HfcP+ Rfe) 

Pfc|fc (Igxg Pfc|fc_i (16x6 K^H/j) -|- 

Lastly, the FSS DKF is considered to be converged when the absolute value of the 
change in the square of the Frobenius norm of the measured diagonal elements of the 
error covariance estimate is below some user-defined threshold, eoKp. This condition is 
shown in Eq. 3.13, where 

(3.12) 

denotes the Frobenius norm of a matrix A: 

I|HfcPfc|fcH^|I p — I|Hfc_iPfc_i|fc_iH^_p|I p < epj^p (3.13) 

While the filter convergence condition in Equation (3.13) is met, the associated bit field 
in the DKF health telemetry is set to "converged," or 1. The filter convergence status is 
utilized in two places in the FSS GNC. First, on test startup, the FSS air bearings and 
guidance subsystem are enabled only after the DKF filter is converged. Secondly, the 
state and error covariance estimates are only passed to the outer-loop for use on the 
next DKF cycle only when the filter is converged. 

3.5.3 Thruster Mapping 

The control subsystem of every FSS includes thruster mapping algorithm. The purpose 
of thruster mapping is to select the appropriate thruster to fire in order to realize the 
desired control input consisting of forces and torques. 

In the development of the thruster mapping algorithm, it is assumed each of the eight 
FSS thrusters produce identical force. The requested control input can be mapped to 
the appropriate thruster from: 

u = Mf (3.14) 

where u = [F^, Fy, r]^ G denotes the control input, f = [/i, / 2 ,..., /s]^ e the 
force applied by each thruster, and M g is the thruster to control input mapping 

matrix. Therefore, the force required by each thruster to realize the control input is 

f = 2M+u (3.15) 
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where M+ denotes the Moore-Penrose pseudoinverse of M. Since the Nf > Nu, M+ 
is given as follows: 

M+= (MM^)“^ (3.16) 

Since the mapping matrix, M, includes both positive and negative assignments, the 
pseudoinverse results in equal scaling along both the positive and negative directions. 
As a result, the factor of two in Equation (3.15) compensates for this fact in order to 
scale the requested force equally across only the applicable thrusters. For instance, if 
a force is requested along the +x direction, thrusters 2 and 3 would be selected, each 
applying half of the requested force. 

Lastly, unlike other actuators which offer a nearly continuous range of control inputs 
across their operating envelope (e.g., reaction wheels, control moment gyroscopes, 
and magnetic actuators), thrusters are discrete and can only take on two states: ON 
or OFF. As a result, difficulty arises where a continuous control signal has to be modu¬ 
lated in a way to actuate the thrusters in order to realize the control signal. A detailed 
description of various thruster modulation techniques is presented in Chapter 4. The 
default modulator, however, is set to be a Sigma-Delta Modulator (SAM). 

3.5.4 Reaction Wheel Controller 

The RW is controlled though a simple speed-mode controller, illustrated in Figure 3.8. 
A speed-mode controller can compensate for the RW friction, and it controls the RW 
such that the correct amount of angular momentum is transferred between the RW and 
FSS. In such a controller, the requested torque is converted, using the RW inertia, to a 
required flywheel acceleration and integrated to obtain a desired flywheel velocity. The 
measured flywheel velocity, ^Measured, is compared with the desired velocity 


Acu — t^desired ^^Measured (3.17) 

and is scaled by a proportional gain in order to generate the required analog voltage to 
drive the RW motor torque. The flywheel velocity is sampled and the RW motor input 
voltage is updated at a 20 Hz rate. To avoid the zero-crossing problem, the RW is spun 
up to about 1100 RPM - half of its maximum rated flywheel velocity - before the air 
pads are enabled and the experiment starts. 


3.6 POSEIDYN Test Bed Characterization 

This section will the detail the system identification of a 4th generation FSS along. 
The characterization includes OS latency, sensor noise, thruster performance, mass, 
moment of inertia, and test bed residual acceleration. 
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^Measured 



Voltage 

Figure 3.8. Reaction wheel speed mode controller. 


3.6.1 OS Latency 

By applying the RT-Preempt patch to the Linux kernel, the maximum latency the OS 
exhibits in response to a stimuli is bounded. To quantify the latency characteristics 
of both the stock (unpatched) Linux kernel and RT-Preempt patched Linux kernel, the 
Cyclictest was run for 50 million iterations [116]. The results are summarized in 
Table 3.3 and illustrated in Figure 3.9. Compared to the stock Linux kernel whose 
maximum observed latency during the Cyclictest was over 13 ms, the PREEMPT-RT 




(a) Stock, unpatched Linux kernel (b) RT-Preempt patched Linux kernel 

Figure 3.9. Comparison of operating system latencies between an un¬ 
patched and patched Linux kernel. 


Table 3.3. Summary of the CyclicTest latency results. 


Kernel Type 

Min Latency 

(/iS) 

Avg. Latency 

(as) 

Max Latency 

(as) 

Stock 

Thread 0 

9 

20 

13462 

Thread 1 

9 

21 

294 

Patched 

Thread 0 

7 

26 

86 

Thread 1 

8 

31 

92 
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patched kernel bounded the OS latencies to within 100 fxs. Furthermore, the maximum 
latency for both the stock and patched kernel was observed to occur on Thread 0. It is 
worthwhile to note, while the average latency of the patched kernel is slightly larger in 
value than the stock kernel, it is of the same order of magnitude and more importantly, 
the maximum latency is bounded to within 100 /xs 

Given the observed average latencies, this corresponds to a latency test duration of 
approximately 17 minutes for the stock Linux kernel and 26 minutes for the patched 
Linux kernel. Assuming an average upper bound on a candidate GNC algorithm test of 
five minutes, the duration of the latency test (for the patched Linux kernel) equates to 
an effective oversampling factor of 5x. Therefore, the TRIDENT-GNC architecture can 
guarantee that the GNC testing code will respond within 0.1 ms — a claim no other test 
bed in literature can make. 


3.6.2 Sensor Noise Characterization 

Vicon Noise Characterization 

It is assumed the noise associated with the measurements provided by the Vicon 
(x,y, and, 6*) are independent since only processed measurements are provided. This 
assumption enables the characterization of each measurement produced by the Vicon 
to be performed independent of the other measurements. To characterize the transla¬ 
tional sensor noise, the FSS is given only a translational velocity across the diagonal 
of the granite monolith. Its position from the Vicon sensor is then recorded directly, and 
the resulting time histories for x and y are independently curve fit to both a first-order 
model, xi(t) and second-order model, X 2 (t), respectively: 

xi{t) = ait + ttQ (3.18a) 

X2{t) = a2t^ + ait + ail (3.18b) 

where { 00 , 01 , 02 } are the model coefficients. Assuming the acceleration coefficient in 
the second-order model, 02 , is sufficiently small, the root mean square error (RMSE) 
of the first-order error is interpreted as the noise parameter for the Vicon x and y mea¬ 
surements respectively. Furthermore, in order to capture any effects of speed on sen¬ 
sor noise, the FSS was given three different speeds across the diagonal of the granite 
monolith. Lastly, the geometric mean was used to compute the overall Vicon noise pa¬ 
rameter for the X and y measurements. A similar method was followed to characterize 
the attitude measurement 6, except, instead of a translational-only velocity, the FSS 
was given only a rotational velocity. 

The resulting Vicon measurement noise values across each trial and the overall geo¬ 
metric means are tabulated in Table 3.4. As tabulated, the Vicon sensor noise for the 
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X measurement is approximately 35% smaller than the y measurement. This is due to 
the arrangement of the Vicon cameras with respect to the granite monolith. The sensor 
noise for the 9 measurement is over 38% larger than that of the y measurement and 
70% larger than the x measurement. This is expected as the attitude measurements 
are derived from the position of each of the markers located atop the FSS, which are 
used in forming the dextral orthonormal basis of the Vicon measurement frame for the 
FSS in the Vicon software. 

Table 3.4. Vicon sensor noise characterization results. 


Trial 

(m) 

ay (m) 

ae (rad) 

1 

0.0113 

0.00889 

0.0356 

2 

0.0118 

0.0130 

0.00255 

3 

0.00199 

0.00662 

0.0269 

Mean 

0.00641 

0.00914 

0.0135 


FOG Noise Characterization 

To estimate the FOG noise parameters, an Allan variance analysis was performed 
following the IEEE Std 952-1997 [117]. The Allan variance method was initially used to 
investigate the frequency stability of precision oscillators in the time domain. However, 
this method has been extended to aid in characterizing random noise processes in a 
wide variety of devices [117], [118] 

To perform the static characterization of the FOG, first, over 13.75 million samples were 
collected at 100 Hz over 27 hours from the FOG, as illustrated in Figure 3.10a. A basic 
analysis of the of static FOG data yielded the mean output, the standard deviation, as 
well as the quantization level. While not relevant to the FOG noise characterization, 
the quantization level is useful for creating the simulated FOG data. It can be seen the 
FOG is able to detect the rotation of the Earth, whose value at the location of the test 
bed is 0.002490 deg Is. 

Next, performing the Allan variance analysis yields the Allan variance plot for the FOG 
illustrated in Figure 3.10b. The relevant noise parameter estimated from the Allan vari¬ 
ance analysis is the angular random walk (ARW). ARW is identified on the Allan vari¬ 
ance plot as having a slope of -1/2. Furthermore, the value for the ARW is obtained 
by fitting a line through this segment on the plot and reading the value at r = 1 [117] 
Given the FOG operates at 100Hz, the equivalent white noise value using the Allan 
variance method was found to be 7.083 x 10“^ deg Is. It is worthwhile to note that the 
white noise value derived from the Allan variance analysis is within 2.3 % of the value 
derived from a statistical analysis of the FOG data. Another useful noise parameter 
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Angular Rate (7s) 


(a) FOG data histogram 

Figure 3.10. FOG data histogram 
and Allan variance plot. 



mean output and la indications 


which can be estimated via Allan variance analysis is the bias stability. This parameter 
would appear on the plot as a segment with a zero slope. Flowever, due to the large 
averaging times required, estimating this parameter via this method is infeasible. Since 
this value is typically smaller than the ARW and given the short duration of the tests, 
this parameter can be neglected. 

A summary of relevant FOG parameters from the two analyses is listed in Table 3.5. 
The white noise value derived from the statistical analysis of the static FOG data to use 
in the DKF was determined to be 7.247 x 10“^ deg Is, or equivalently, 1.265 x 10“^ 
rad Is. 


Table 3.5. Summary of FOG noise analysis. 


Parameter 

Value 

Mean Output 

1.813 X 10“^ deg /s 

White Noise 
(Statistical Analysis) 

7.247 X 10-3 deg /s 

White Noise 
(Allan variance) 

7.083 X 10-3 deg /s 

Ouantization Level 

1 X 10-® deg /s 
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3.6.3 Thruster characterization 

Characterizing the performance of the thrusters is essential across several aspects of 
the POSEIDYN test bed. This includes accurate modeling in the numerical simulator, 
control input in the DKF, guidance parameter tuning, and moment of inertia (MOI) es¬ 
timation. A single thruster was characterized in [100]; here, a characterization of the 
thruster performance when integrated into the FSS is presented. Given the architecture 
of the system, the force generated by the thruster is highly influenced by the pressure 
regulator setting and the number of values open at the same time. In an effort to char¬ 
acterize the nominal performance of the thrusters, two thrusters on the same face are 
fired for several seconds and the maximum force generated by the thrusters is mea¬ 
sured by a scale with a resolution of 9.81 mN. Additionally, it is assumed the measured 
force has equal contributions from the two thrusters. Prior to firing the thrusters, the 
value reported by the pressure regulator, with a one PSI (0.0689 bar) resolution, is 
recorded such that a pressure versus thruster force mapping can be made. This pro¬ 
cess is repeated across a range of operating pressures and several times to ensure 
repeatability. The results are illustrated in Figure 3.11. As illustrated, several regula¬ 
tor pressures took on similar thrust values due to the resolution of both the scale and 
pressure regulator. 



Figure 3.11. FSS thruster characterization data with first-order fit and es¬ 
timated 3(j bounds. 


Since no other data, such as temperature, is available, a linear fit is made between the 
regulator pressure and thrust generated. This linear model can then be used in the 
MOI estimation as well as adjusting the thrust parameter for the DKF and simulator. It 
is important to note, the onboard GNC software does not receive pressure information 


44 








from the regulator which precludes the use of the linear model onboard for estimating 
thrust. 


3.6.4 Mass and Moment of Inertia 

Estimating the physical properties (i.e., mass and MOI) is critical to not only under¬ 
standing the system dynamics, but also to creating a realistic simulator and developing 
and tuning the guidance algorithms. The system masses of a FSS are tabulated in 
Table 3.6. 

To estimate the FSS MOI, a torque was first applied to the stationary vehicle in the -9 
direction by the thrusters, followed by a coast period; lastly, an opposite torque was 
applied in the +9 direction by the thrusters to spin the vehicle down. It is worthwhile 
to note, the value indicated by the thruster pressure regulator was recorded prior to the 
test in order to use the linear model (from Figure 3.11) to estimate both the thrust output 
and torque of each thruster. The resulting attitude time history was recorded directly 
from the Vicon sensor. Next, the portions of the attitude time history under thrusting 
were fit to the second-order model 

9(t') = (22^^ “h Q-it Oq (3.19) 

where 02 = l/2(r/J^;j), from which the MOI (J^z) is estimated given the estimated 
torque r. This test was performed five times in order to gather 10 attitude time histories 
under thrust, refueling the FSS each time prior to running the test. The resulting MOI 
estimate is tabulated in Table 3.6. 

Table 3.6. Summary of relevant FSS physical properties. 


Parameter 

Value 

Mass, wet 

9.882 kg ±0.001 kg 

Mass, dry 

9.465 kg ±0.001 kg 

Dimensions 

0.27x0.27x0.52 m 

Estimated MOI 

0.253 kg-m^ 

MOI Error (la) 

0.0115 kg-m2 


3.6.5 Residual Acceleration 

From an experimental standpoint, it is important to develop a bounded estimate of the 
residual accelerations acting on the FSS during a test campaign. In addition to any 
component of gravity in the plane of motion and residual friction, these unmodeled ac¬ 
celerations are also due to lateral air bearing leakage due to the air bearing not exactly 
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parallel to the test bed surface and the air flow coming out of the air bearing being not 
exactly symmetric [92], In an attempt to quantify these residual accelerations, a FSS 
was positioned at the nine locations illustrated in Figure 3.12 on the granite monolith 
via an onboard Linear Quadratic Regulator (LQR) controller. Upon arriving at a pre¬ 
defined grid point, the air pads were closed for 30 seconds, allowing the FSS and the 
DKF to stabilize, before opening again for 30 seconds to perform the test at the loca¬ 
tion. The resulting acceleration imposed on the FSS during the 30 second interval was 
then derived from on-board state estimate information. The average (scaled) acceler¬ 
ations at each location over four tests is illustrated in Figure 3.12. The overall average 
translational and rotational residual acceleration over the nine locations was found to 
be 1.871 X 10“'^m /s^ (19.1/iG) and 7.56 x 10“^ deg /s^, respectively. 



X-Position (m) 

Figure 3.12. Average residual acceleration at nine locations on the granite 
monolith. 


3.7 Case Study 

To demonstrate the capabilities of the POSEIDYN test bed, an exemplary case study 
is performed where the FSS enters into and maintains a circular path with a 0.50 m 
radius and a tangential velocity for four revolutions, at a constant angular velocity of 3 
deg /s about the center, using an APF-based guidance law to control the FSS [19], [20], 
[28], [32]. The APF-based guidance law was chosen as it is a real-time guidance and 
control method. 
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3.7.1 Guidance and Control Development 

Assuming some initial configuration of the FSS that is outside of the desired circular 
reference path, the guidance will initially target a point, p/, that is tangential to the 
current position of the FSS, p, and the circular reference trajectory in the direction of 
the desired rotation around the center: 


P/ 


rref sign 




+ Pref 


(3.20) 


where net is the radius of the circular reference trajectory; ^/>ref is the reference angular 
velocity about the center of the circular trajectory; pref = [xc-,y^ is the center of the 
reference trajectory; /3t = sin“^ {{r^Q^/D)) is the angle between the FSS and a point 
tangent to the reference trajectory and D is the distance from the FSS to the center 
of the reference trajectory; and 0 = tan“^ ((?/c - y)/{xc - x)) is the angle of the FSS 
relative to the center of the reference trajectory in the inertial frame. During this seg¬ 
ment, only the target, or desired, position, p/, is passed to the APF-based control law. 
Furthermore, the FSS is commanded to point the +x face towards the tangent intercept 
point. It is worthwhile to note, attitude control is performed via thruster actuation. 


Once the FSS is within a user-defined distance, to the targeted point p/, the guidance 
then targets the circular reference trajectory and sets the initial angle 0o, the angle 
of the FSS relative to the center of the reference trajectory with respect to the inertial 
x-direction at this point in time. Subsequent reference positions and velocities are 
determined by propagating the initial reference angle at the user-specified rate, 0ref: 


P/ 

P/ 


Pref + r-ref 


COS (0) 
sin (0) 


r’ ref ^ref 


— sin (0) 
COS (0) 


(3.21) 

(3.22) 


During the circumnavigation portion of the maneuver, the +x face of the FSS is pointed 
inwards towards the fixed reference point 


^'ref 


tan ^ 


/ yc-y 

\Xc - X 


(3.23) 


An APF control law is used to generate the desired commands to follow the reference 
trajectory and reference attitude, x/. To generate the commands, a potential field is 
setup utilizing only one quadratic attractive potential located at the reference position 
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and attitude 


(3.24) 


- X/)^Qa(x - X/) 

where e is a symmetric, positive-definite gain matrix, and Qa e is 
the symmetric, attractive potential shaping matrix. It is worthwhile to note, if additional 
obstacles or keep-out zones were present, repulsive potentials could be placed at these 
locations to ensure obstacle avoidance and keep-out zone enforcement. The resulting 
control law is given as the gradient of the attractive potential (and repulsive potential, if 
used): 

U = (Qa(x - Xref) + (x - Xrgf) (3.25) 

where B e is the lower half of the continuous-time input gain matrix. 


3.7.2 Case Study Setup 

A list of the relevant test parameters are reported in Table 3.7. It is worthwhile to note, 
prior to running the test on the POSEIDYN test bed, the typical development workflow 
included preliminary testing using the development numerical simulator to tune and 
debug the guidance algorithm. 


Table 3.7. Case study test parameters. 


Guidance 

Parameter 

Value 

Navigation 

Parameter 

Value 

Control 

Parameter 

Value 

Sample 

Time 

0.10 s 

Sample 

Time 

0.02 s 

Sample 

Time 

0.01 s 

{Xrei, yref) 

(2, 2) m 

4 

0.001 

SDM Sample 
Time 

0.10 s 

?’ref 

0.5 m 


6.4123 mm 

Thruster 

Force 

0.115 N 

i’rei 

3 deg /s 

CJy 

9.1363 mm 

Min. 

Firing Time 

0.06 s 

No. of 

Revs 

4 

ere 

1.3451 mrad 

Min. 

Impulse Bit 

0.0115 N s 

Qa 

diag(0.05,0.05,0.05) 

^e 

0.12649 mrad /s 

— 

— 

Ka 

diag(1,1,0.1) 

c 

12.116 

— 

— 

— 

— 

.2 

^DKF 

5 X 10“® 

— 

— 
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3.7.3 Case Study Results and Discussion 

As illustrated in Figure 3.13 and Figure 3.14, the FSS starts in the lower-right hand 
corner of the granite monolith. At the start of the test, the guidance subsystem, tar¬ 
gets a point tangent to the circular reference trajectory in the direction of the reference 
angular velocity. It can be seen that during the tangent point targeting segment the 
FSS attitude only achieves a near-zero attitude error as it nears the reference trajec¬ 
tory transition (marked by a red asterisk, *) due to a small attractive potential weighting 
on the attitude error. Once the FSS gets sufficiently close to the circular reference tra¬ 
jectory, the reference trajectory switches to the and the FSS is seen to track the desired 
circular trajectory fairly well. As illustrated in Figure 3.15, the maximum distance from 
the circular reference trajectory was less than approximately 0.04 m. The maximum 
difference between the FSS speed and the reference error speed was no larger than 
approximately 0.003 m/s, except for a brief period of time while the FSS was getting 
established on the reference trajectory. Furthermore, a relatively constant attitude error 
throughout the circumnavigation can be seen in Figure 3.13 and Figure 3.15c. This 
is attributed to the small attractive potential weighting on the attitude error used in the 
case study. The time, position, speed, and attitude errors across each revolution and 
the entire maneuver are summarized in Table 3.8. Lastly, as anticipated, frequent ac¬ 
tuated control inputs occur throughout the maneuver to maintain the circular reference 
trajectory as illustrated in Figure 3.16. 



X-Position (m) 


Figure 3.13. Overhead view of the FSS path with four equally-spaced atti¬ 
tude snapshots in time across the tangential intercept guidance segment 
and a single revolution about the fixed point. (Traj: Trajectory) 
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Position (ra) 




(a) Position (b) Velocity 

Figure 3.14. Position and velocity of the FSS during the entire maneuver. 


Table 3.8. Summary of the circumnavigation time, position, speed, and 
attitude errors. 



1 

Revolution 

2 3 

4 

Overall 

Rev. Time, s 

120.00 

120.00 

120.00 

119.91 

119.98 

Mean Pos Error, mm 

9.18 

6.56 

6.64 

6.64 

7.25 

Mean Speed Error, mm/s 

1.18 

1.01 

1.01 

0.99 

1.04 

Mean Att. Error, deg 

35.07 

26.64 

27.62 

27.01 

29.08 

Mean Att. Rate Error, deg /s 

0.52 

0.34 

0.35 

0.30 

0.38 


Table 3.9. Mean estimated error standard deviation over the circumnavi¬ 
gation maneuver. 


Revolution 



1 

2 

3 

4 

Overall 

X, mm 

2.26 

2.63 

2.64 

2.65 

2.63 

Va;, mm/s 

4.19 

4.20 

4.21 

4.21 

4.21 

y, mm 

3.39 

3.41 

3.43 

3.44 

3.42 

Vy, mm/s 

4.58 

4.59 

4.60 

4.60 

4.59 

0, deg 

0.037 

0.038 

0.036 

0.036 

0.037 

9, deg Is 

0.016 

0.016 

0.015 

0.014 

0.015 
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o to *>. 




w E 0 
^ -4 


0 0.5 1 1.5 2 2.5 3 3.5 4 

Revolution Number 

"O 0.5 1 1.5 2 2.5 3 3.5 4 

Revolution Number 

4 I- ^^^^^^^-1 

0 0.5 1 1.5 2 2.5 3 3.5 4 

Revolution Number 


(a) Position Error 


(b) Velocity Error 




1.5 2 2.5 

Revolution Number 


(d) Attitude Rate Error 

Figure 3.15. Position and velocity errors of the FSS during circumnaviga¬ 
tion. 
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(a) F, 


(b) fy 



(c) Torque (thruster actuation only) 


Figure 3.16. Commanded and actuated control inputs (in the inertial 
frame) throughout the circumnavigation maneuver. 
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Attitude Covariance (°) 



(c) y Covariance 


(d) Vy Covariance 




Revolution Number 


(f) Attitude Rate Covariance 


Figure 3.17. Covariance time histories throughout the circumnavigation 
maneuver. 
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Figure 3.18. Detailed illustration of the Attitude Rate Covariance recovery. 


To gauge the performance of the FSS DKF, both the navigation health telemetry and 
the estimated error covariances must be examined. The estimated error covariances 
throughout the duration of the maneuver is illustrated in Figure 3.17 and is summarized 
in Table 3.9. As a result of the frequent thruster firings required to track the desired tra¬ 
jectory, as illustrated in Figure 3.16, this case study tests the performance of the DKF. 
Specifically, the DKF is given minimal time to filter the noise introduced by the thruster 
firings. Additionally, the effects of increased process noise on the attitude rate covari¬ 
ance is illustrated in Figure 3.18. Recall, whenever a thruster is actuated, the navigation 
filter automatically inflates the process noise to account for uncertainties in the thruster 
model. Despite the process noise associated with the attitude rate state being inflated 
by a factor of 500, the attitude rate covariance is able to recover quickly (within 1 to 2 
DKF cycles) due to the small sensor noise associated with the FOG. Lastly, the small, 
periodic, increases in the attitude rate covariance in between thruster firings were due 
to unavailable FOG sensor data, as indicated by the DKF health telemetry. 


3.8 Summary 

A detailed description of the Naval Postgraduate POSEIDYN test bed has been pro¬ 
vided. The quasi-frictionless and low residual acceleration dynamic behavior of the test 
vehicles operating on the POSEIDYN test bed are limited to two translational and one 
rotational degree-of-freedom motion. Additionally, multiple vehicles can be operated 
simultaneously on the POSEIDYN test bed. This capability can be used to conduct re- 
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search on coordinated control of spacecraft teams. The on-board actuators, composed 
of eight thrusters and/or a reaction wheel, are equivalent to the actuators found in orbital 
spacecraft and further enhance the dynamic equivalence. With an available simulator 
and extensive software development tools, guidance and control algorithms for real¬ 
time execution on-board the test vehicles can be quickly and easily developed. Used 
as a last stage of on-the-ground validation prior to on-orbit deployment or simply as 
a more realistic development environment for novel guidance and control approaches, 
this state-of-the-art dynamic hardware-in-the-loop test bed will continue to be fruitful for 
the advancement of spacecraft proximity operations research. 
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CHAPTER 4: 

Spacecraft Thruster Modulation Techniques 


Delta Sigma Converters or Sigma Delta Converters? Mankind 
does not seem to agree on one notation. Both notations are used 
equally often when you search via Google. I decided to stay with 
that guy who told he Is living In the Mississippi Delta, so deltas 
mean something to him - and for him only the Sigma River may 
have a Sigma Delta...good point. Later I found out that the original 
name “Delta Sigma” was coined by the Inventors Inose and 
Yasuda and "Sigma Delta Is actually not correct. I was lucky... 

— U. Beis 

Quoted by Hebertt SIra-Ramfrez 
Sliding Mode Control [119] 


A journal version of the work presented in this chapter has been accepted for publication 
in the AIAA Journal of Guidance, Control, and Dynamics [70]. 


4.1 Overview of Thruster Modulation Techniques 

Spacecraft reaction control (RC) thrusters provide a means for attitude control as well 
as orbital maneuvering. When compared to other actuators, offering a nearly contin¬ 
uous range of control inputs across their operating envelope (e.g., reaction wheels, 
control moment gyroscopes), RC thrusters take on only two discrete states: ON or 
OFF. To use the RC thrusters, the often continuous control signal has to be modulated. 

One of the simplest thruster control methods is the Schmitt Trigger, commonly de¬ 
scribed as a relay with hysteresis and a deadband [82]. This control technique, which 
is not technically a pulse modulator, defines its minimum pulse-width as a function 
of the changing spacecraft’s inertia [82]. The Pulse-Width Pulse-Frequency Modula¬ 
tor (PWPFM) and Derived-Rate modulator extend the Schmitt Trigger with the addi¬ 
tion of a first-order lag filter in either the feed-forward or feedback paths respectively. 
Another alternative is the Pulse-Width Modulator (PWM), resembling a PWPFM in be¬ 
havior, but much simpler in its construction. The static characteristics of the PWM, 
PWPFM, and Derived-Rate modulator are independent of the spacecraft inertia [82], 
which is an advantage over the Schmitt Trigger. The PWM, PWPFM, and Derived-Rate 
modulator have been widely used for spacecraft thruster modulation [120H122]. 
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Here, the use of the Sigma-Delta Modulator (SAM)^is investigated and proposed as an 
alternative modulation technique for spacecraft thruster modulation. Originally applied 
to analog-to-digital converters by Inose and Yasuda, the SAM has been employed 
in numerous applications, ranging from communications, analog-to-digital converters, 
high-quality audio playback, to high-resolution instrumentation [123H125]. Ciarcia et 
al. briefly considered the use of a SAM for the RC thruster modulation for the trans¬ 
lational control of a spacecraft simulator — however, due to its estimated control cost 
(in an ideal simulation environment), a different actuation strategy was chosen for their 
application [69]. To the author’s best knowledge, the use of a SAM for RC thruster 
modulation has not been previously used. 

Compared to the preceding pulse-control modulators (PCMs), the SAM utilizes a coarse 
(typically 1-bit) quantizer in combination with feedback to realize a “high resolution 
conversion of a relatively low bandwidth signal” [126]. Furthermore, this modulation 
technique provides user-defined noise shaping as well as small-signal actuation [124]- 
[126]. When compared to other PCMs, the SAM is capable of realizing smaller input 
signals over smaller sampling periods. The SAM can also be linearized through multi¬ 
ple methods, enabling the application of traditional continuous-time design and analysis 
techniques [119], [123], [125], [126]. 

Given its unique features and apparent advantages over traditional PCMs, the use 
of a a first-order SAM as a thruster modulation technique is examined in detail and 
experimentally compared to a PWM. A PWM is chosen for the comparison for its 
simplicity and wide adoption on spacecraft RC control [127], [128]. 

Linearized models of each modulator are presented first. Next, these models, in combi¬ 
nation with a simple PD control law, are used to develop candidate closed-loop models. 
An experimental campaign is then performed using the NPS POSEIDYN test bed, com¬ 
paring the PWM and SAM in a relevant environment with hardware-derived noise and 
uncertainties [31], [33]. Finally, the results from these experiments, and from various 
numerical simulations, are analyzed to draw conclusions regarding the applicability of 
the candidate models and of the SAM’s merit as a thruster modulation technique. 

^ Despite the anecdotal epigraph at the beginning of the Chapter, the term “Sigma-Delta” will be used 
throughout this Chapter. According to Gray, the alias “Sigma-Delta” is attributed to “[J.C.] Candy and 
his colleagues” [123]. A common argument in-favor of the use of the alias is that it is a more accurate 
description of the modulator, as the summing (S) is performed over the differences (A) [123], [124]. 
Additionally, Gray further argues that “the system does not really form a difference of successive samples 
[such as the A modulator] of the input signal as suggested by the AE name” (p. 60) [123]. 
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4.2 Pulse Modulator Development 

The construction of a first-order SAM and PWM and their linearized models is pre¬ 
sented in this section. These models are later used to derive candidate models of 
an attitude control system with a modulator-in-the-loop (MITL). These functional de¬ 
scriptions can either be viewed from either a hardware (e.g., analog-circuit) or digital 
implementation (e.g., software) perspective, but emphasis is placed on digital imple¬ 
mentation, as this is the typical implementation medium. 

4.2.1 Sigma-Delta Modulator 

Overview 

The topology of a first-order SAM, illustrated in Figure 4.1 a, takes the form of negative 
feedback loop. Along the feed-forward path, a low-pass filter, in this case a single 
integrator, is in cascade with a 1-bit quantizer; the feedback path is composed of a 
1-bit digital-to-analog converter [124]-[126]. Without loss of generality, the input to the 
SAM is a time-varying signal, v(t). In this context, v(t) represents the commanded 
force to be generated by a particular thruster. The input signal is sampled via a zero- 
order hold (ZOH) mechanism at a frequency /eam = 1/7sam [129]. The output of 
the modulator, y(t), is a digital pulse train and is converted back to an analog signal 
by the 1-bit digital-to-analog converter along the feedback path. In this application, 
the feedback signal can consist of either a model-driven or sensor-driven response of 
the thruster. As such, the sampling period of the SAM should be similar to the impulse 
period of the thruster. The error signal, e(t), is generated by subtracting the output from 
the input. It is then multiplied by a strictly positive gain k eR+ before being integrated 
over the sampling period T^am to form the intermediate signal w(t) = k f e(t) dt. Lastly, 
the intermediate signal, w(t), is modulated via a 1-bit quantizer to form the output 
signal y(t), which is then fed-back. A pulse is triggered when w(t) reaches a desired 
threshold [123], [125]. In this particular application, the quantizer threshold, e, is lower- 
bounded by the thruster minimum impulse bit. 


Linear Model Development 

First, it is assumed the modulator sampling frequency, is sufficiently fast such that the 
ZOH can be neglected, allowing the SAM to be viewed as a continuous-time device. 
The only remaining nonlinear element in the SAM is the quantizer, which cannot be 
directly converted into an equivalent linear system. Multiple approaches can be taken 
which will linearize this modulator under various assumptions. The first method, which 
has been developed and used throughout the 20th century, assumes the quantizer er¬ 
ror can be adequately described by an additive white noise source, as illustrated in 
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Figure 4.1. First-order EAM topology and linearized block diagram. 


Figure 4.1b. This assumption, commonly known as the additive white-noise approxi¬ 
mation assumes the quantization error is uniformly distributed over the two quantization 
levels provided by the 1 -bit quantizer [123], [125]. In order to satisfy this assumption, it 
is necessary that the magnitude of the input to the modulator, |f(t)|, remains bounded 
by the output of the modulator, Vmax [123], [125]. This results in a two-input, one-output 
closed-loop linear system given as follows: 


Y(s) 


k 

s k 


VW + 


s 

s + fc 


N(s) 


(4.1) 


Note, the high-pass noise-shaping property of the EAM is highlighted by the noise 
transfer function (NTF) given by 


The corresponding magnitude response of the noise transfer function is 

00 

HntfOi^) = 



(4.2) 


(4.3) 


which is that of a high-pass filter with a cutoff frequency ot oo = k since IHntfO'^)! = 
1/a/2. Resultantly, the EAM effectively “pushes” the noise to higher frequencies while 
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simultaneously increasing the signal-to-noise ratio for frequecies below the cutoff fre¬ 
quency. 

Another method used to linearize the EAM views the modulator as a sliding mode con¬ 
troller where the quantizer is represented as an ideal relay [119]. In this formulation, the 
quantizer is represented as a relay with an infinite switching bandwidth - thus convert¬ 
ing the modulator to a continuous-time device. By viewing the SAM as a sliding mode 
controller, the average output of the modulator is equivalent to its input. However, since 
the output of this model is a switched-output, a low-pass filter is required to demodulate 
the signal into a continuous waveform [119]. 

4.2.2 Pulse-Width Modulator 

Overview 

In contrast to the SAM, the PWM is an open-loop modulator, whose topology is illus¬ 
trated in Figure 4.2a, that generates a series of varying-width pulses that represent the 
input signal [130]. Similar to the SAM, the input to the PWM is a time-varying signal, 
v{t), representing the desired thrust level, which is sampled, via a ZOH, at the PWM 
frequency, /pwm = [129], [130]. Next, the input signal is normalized by the 

maximum output of the PWM, in this case, the thruster output, V^ax, forming the de¬ 
sired pulse width, ^(t) = v{t)/Vmax- To realize the pulse-width, a normalized repeating 
reference waveform r(t), with period Tr, is sampled at an integer n G M+ frequency 
greater than the PWM frequency, /r = n/pwm [130]. Note, typical reference waveforms 
used include a Sawtooth, Inverted Sawtooth, and Triangle waveform [130]. During a 
given PWM period, the output of the PWM, y{t), is a digital pulse whose width which 
corresponds to the amount of time when ^{t) > r{t). 




Modulator 

Output 


Signal 


Comparator 


(a) PWM Topology 



(b) Linearized PWM 

Figure 4.2. PWM topology and linearized block diagram. 
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Compared to the SAM which has only two discrete output states (ON or OFF), the 
number of discrete output levels of the PWM is controlled by the integer ratio 


pwm 


T 

pwm 


(4.4) 


Note, the sampling period of the normalized repeating reference waveform, deter¬ 
mines the pulse width resolution, which establishes a lower bound on the deadband 
associated with this modulation technique [82], The upper bound, which can be dif¬ 
ferent from the minimum pulse width, is given by the minimum impulse of the thruster. 
Note, the lower bound is affected by the implementation of the PWM (e.g., Tr), whereas 
the upper bound is dependent upon the actuator hardware (e.g., thruster solenoid re¬ 
sponse time). 


Linear Model Development 

The only nonlinearity present in the PWM is the comparison operation to determine the 
pulse-width. Since the fr » /pwm, the reference signal ZOH can be neglected. How¬ 
ever, since the PWM frequency is much smaller than the reference signal frequency, 
but is still larger than twice the largest frequency of concern in the input signal, the 
PWM ZOH may not be neglected. Viewing the comparison operation as a quantization 
process, the same approximation can be made regarding the linearity of the quantiza¬ 
tion error. Assuming the input v(t) is bounded similarly to the SAM, the comparison 
operation can be viewed as a quantization process and the same approximation can 
be made regarding the linearity of the quantization error. The resulting linearized PWM 
model, as illustrated Figure 4.2a, becomes a ZOH with a sampling period of Tpwm, 
whose transfer function is 


Y(i>) 


g—Tpvvm^ 


s 


V(s) + N(s) 


(4.5) 


4.3 Closed-Loop Spacecraft ACS with a Modulator-in- 
the-Loop 

A pertinent application of spacecraft thruster actuation is for attitude control. Given 
the notional spacecraft attitude control system illustrated in Figure 4.3, the commanded 
actuation from the attitude controller is passed to the steering logic which selects the 
appropriate thrusters. The resulting desired thruster output signal is then modulated 
,which, in-turn, drives the actuator. To gain insight on the effect of each modulation 
method, a closed-loop model describing the behavior of an attitude control system with 
a MITL will be developed. 
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Ground 

Commands 



Figure 4.3. Notional spacecraft attitude control system block diagram. 


The attitude dynamics of a spacecraft follow Euler’s equation 

T = JcJ (4.6) 

where r G is the net torque acting on the spacecraft, J G is the inertia matrix 
with respect to the COM in the spacecraft body Cartesian coordinate system; cu g is 
the angular velocity of the spacecraft in the body frame with respect to the inertial frame 
expressed in the spacecraft body Cartesian coordinate system, and is the matrical 
representation of the cross product operator [82]. When the inertia is measured in the 
Principal Axis Cartesian coordinate system, the inertia matrix becomes diagonal. For 
small angular velocities, the gyroscopic term becomes negligible and Equation (4.6) 
becomes, 

T = Jdj (4.7) 

This implies the motion about each axis is decoupled and can be controlled indepen¬ 
dently. Furthermore, assuming no external disturbances, r becomes the control torque. 
For simplicity, a single-axis, double integrator model representing the spacecraft atti¬ 
tude dynamics is considered here. 

Combing the transfer functions for each modulator and the attitude dynamics yields. 


Gpwm('S) — 

_ g—-ipwmS 

(4.8a) 


Gsam(s) = 

1 

(4.8b) 

J(s 3 + s2) 


Assuming a PD control law of the form 

C(s) = kdS + kp (4.9) 
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where kp and kd are the proportional and derivative gains respectively, the candidate 
closed-loop model are 


Tpwm ("5) 
Tsam(s) 


kp + - 1 ) 

(Js^ + kdS + kp) — kdS — kp 

{kd^ “h kp) 

Js^ + + kdS + kp 


The validity of these models will be assessed in Section 4.4.3. 


(4.10a) 

(4.10b) 


4.4 Experimental Comparison 

This section will first discuss the experimental test setup and selection of pertinent 
experimental parameters. The results of the experimental campaign are then presented 
and discussed. Lastly, the validity of the linear candidate models are assessed. 

4.4.1 Experimental Setup 

A first-order EAM was experimentally tested and compared to a PWM to demonstrate 
its performance and efficacy for spacecraft thruster control. A single-axis, 90 deg, 
rest-to-rest, slew maneuver was performed using an underdamped PD controller to 
demonstrate both the signal and steady-state tracking abilities of each modulator. The 
following metrics were leveraged in comparing the two modulators: thruster ON-time, 
steady-state attitude error, steady-state attitude rate error, and number of thruster actu¬ 
ations. Relevant experiment and spacecraft simulator system parameters are summa¬ 
rized in Table 4.1. Each experiment was repeated 10 times on the NPS POSEIDYN test 
bed [31], [33]. The setup for the experimental comparison is illustrated in Figure 4.4. 



Figure 4.4. Test setup for the comparison between EAM and PWM. 
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Before continuing, it is worthwhile to comment on the selection of a few of the experi¬ 
mental parameters. The sampling period of the SAM was chosen to be slightly larger 
than the minimum solenoid response time as to ensure adequate flow development 
from the thruster, while the quantization level was chosen to be the minimum impulse 
bit for each thruster. The PWM minimum impulse time - the minimum amount of time for 
which the thruster must actuate - was chosen to be equivalent to the thruster solenoid 
response time, while the Pulse-Width resolution, T^, was chosen to reduce the dead¬ 
band associated with the attitude control [82]. Using Equation (4.4), the PWM sampling 
period of 1 s was chosen to allow for 100 discrete output levels. 

Table 4.1. Summary of relevant experiment and FSS parameters. 


Parameter 

Value 

Experiment Parameters 

Experiment Duration 

120 s 

Number of Trials 

10 

Reference Angle, Of 

90 deg 

Reference Attitude Rate 

0 deg /s 

Attitude Controller Sample Rate 

10 Hz 

Proportional Gain, kp 

0.10 

Derivative Gain, kd 

0.20 

FSS Parameters 

Mass 

9.946 kg 

Est. Inertia 

0.253 kg-m^ 

Est. Thruster Force 

0.117 N 

Est. Total Torque 

0.0437 N m 

Thruster Solenoid Response Time 

0.06 s 

SAM Parameters 

SAM Sampling Period, Team 

0.1 s 

SAM Error Gain, k 

1 

SAM Quantizer Threshold, e 

0.0117 N s 

PWM Parameters 

PWM Minimum Impulse Time 

0.06 s 

PWM Sampling Period, Tpwm 

1.0 s 

Pulse-Width Resolution, 

0.01 s 

Number of Discrete Output Levels 

100 
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4.4.2 Experimental Results & Discussion 

The results presented in Table 4.2 and Figure 4.6 are averaged over 10 trials. For the 
experimental parameters listed in Table 4.1, the SAM was found to require approxi¬ 
mately, on average, an additional 3 s of thruster ON-time when compared to the PWM 
but achieved a significantly smaller attitude error. This difference in thruster ON-time is 
attributed to the feedback loop in the SAM, resulting in shorter period limit cycle, when 
compared to the PWM. After reaching a steady-state (at approximately 25-35 s), the 
SAM was observed to command, on average, 22.5 thruster pulses (approximately one 
firing every 3.6 s) achieving an average attitude error of 0.32 deg. The PWM, however, 
was observed to actuate the thrusters, on average, 6.50 times throughout a similar 
period (approximately once every 12.3 s) achieving an average attitude error of 2.23 
deg. Both modulation methods produced a nearly zero average steady-state attitude 
rate error; the standard deviation for the PWM was observed to be about half of that of 
the SAM. As illustrated by Figure 4.5a, both modulators produce a similar (average) 
phase plane trajectory. In contrast to the smoother trajectory produced by the SAM, 
the trajectory produced by the PWM has distinct and pronounced coast and actuation 
segments. 


Table 4.2. Comparison of the 10-run averaged experimental metrics (with 
Icr error) between the SAM & PWM. 



SAM 

PWM 

Thruster 

ON-Time, s 

24.20 ± 0.47 

21.13 ± 1.13 

Steady-State 

Attitude Error, deg 

0.32 ± 1.59 

2.23 ± 2.60 

Steady-State 

Attitude Rate Error, deg Is 

0.00 ± 0.73 

0.01 ± 0.36 

Number of 

Thruster Actuations 

100.50 ± 34.85 

49.17 ± 17.99 


Several important conclusions regarding the performance of the two modulators can 
be drawn by comparing the thruster ON-time history against time and the attitude and 
attitude rate errors, illustrated in Figure 4.5b and 4.5c and 4.5d. First, it becomes ap¬ 
parent that the SAM leads the response of the PWM and resultantly settles to within 5 
deg of the reference attitude 9f approximately 10 s faster. As illustrated by the orange 
asterisks (*) in Figure 4.5c, the SAM achieves steady-state using approximately 3 s 
less thruster ON-time! However, due to the short period of the limit cycle induced by 
the feedback in the SAM, the thruster ON-Time was observed to linearly increase at a 
faster rate than the PWM. Given a set of mission requirements, the SAM can be cou¬ 
pled with a deadzone in order to increase the period and size of the limit cycles, thereby 
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Thruster ON-Time (s) Attitude Rate Eiror ( Vs) 


reducing the fuel consumption. Alternatively, one can utilize the SAM throughout the 
transient period and switch to a lower-bandwidth modulator, such as the PWM, when in 
the vicinity of the reference state. 




(a) Phase Plane Comparison 


Time (s) 

(b) Thruster ON-Time Time History 



(c) Attitude Error vs. Thruster ON-Time 



(d) Attitude Rate Error vs. Thruster ON-Time 


Figure 4.5. Direct comparison of the SAM and PWM trajectories averaged 
over 10 runs. 
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Attitude Rate Error (7s) Torque (N-m) Attitude Error (“) 




(a) SAM: Attitude Error Response (b) PWM: Attitude Error Response 



(c) SAM: Control Time History 



(d) PWM: Control Time History 



Attitude Error (°) 

(e) SAM: Phase Plane Response 
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(f) PWM: Phase Plane Response 


Figure 4.6. Experimental results averaged over 10 runs comparing the 
response of the SAM to PWM for a 90 deg slew. 
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4.4.3 Linear Model Validation 

The candidate models developed which describe the closed-loop behavior of a space¬ 
craft attitude control system with a MITL are described by Equations (4.1 Oa) and (4.1 Ob). 
To validate the candidate models, a comparison was made of their (numerically simu¬ 
lated) attitude error time histories against both the (averaged) experimental response 
and a (numerically) simulated nonlinear response of FSS [31], [33]. The idealized (lin¬ 
ear) response with no MITL is presented for reference purposes. When implementing 
the numerical simulations, the control input was saturated at the maximum total torque 
listed in Table 4.1 for all systems. The resulting comparison is illustrated in Figure 4.7. 
Note, only the first 60 seconds are shown in order to increase the clarity of each re¬ 
sponse. Furthermore, the point of peak overshoot is indicated in Figure 4.7 by a dashed 
line (—) while the settling time is noted by a dashed-dot line (- ■ -). Transient and 
settling characteristics of each system are tabulated in Table 4.3. 


Table 4.3. Comparison of transient and settling characteristics. 




Ideal 

Linearized 

Nonlinear 

Experimental 



Response 

Response 

Response 

Response 

Percent Overshoot, % 

SAM 

7.12 

38.94 

45.15 

42.44 

PWM 

7.12 

13.71 

44.89 

42.20 

Peak Time (Tp), s 

SAM 

7.43 

8.15 

14.90 

12.79 

PWM 

7.43 

7.49 

15.10 

14.19 

Settling Time (T*), s 

SAM 

8.67 

18.25 

38.29 

24.70 

PWM 

8.67 

9.32 

39.99 

34.70 

Tp/Ts 

SAM 

0.86 

0.45 

0.39 

0.52 

PWM 

0.86 

0.80 

0.38 

0.41 


While the (simulated) attitude error response of the candidate SAM model did not 
exactly match that of the experimental or simulated nonlinear responses, it did exhibit 
similar characteristics. First, the candidate model was observed to have a percent 
overshoot of 38.94%, corresponding to a 8.6% difference to the experimental response 
and a 14.8% difference to the simulated nonlinear response. Furthermore, comparing 
the ratio of peak time, Tp, to settling time, T^, the candidate model was found to obtain 
a value (approximately) half-way in-between the experimental and simulated nonlinear 
responses, corresponding to a 14.3% difference to the simulated nonlinear response 
and a 14.4% difference to the experimental response. This has the implication of having 
similar damping characteristics to the two comparison responses. Given the similarities 
in the transient and settling characteristics of the candidate SAM model compared 
to the experimental and simulated nonlinear responses, this model is considered to 
adequately describe the physical system. 
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The PWM candidate model, on the other hand, was not found to be in agreement with 
either the experimental of simulated nonlinear responses. This is apparent due to a 
significant mis-match between the transient and settling characteristics of the candidate 
model and the two comparison responses. 



(a) SAM 



(b) PWM 

Figure 4.7. Comparison of the linearized model response to the simulated 
and experimental responses. 
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4.5 Summary 

Experimentally comparing a Sigma-Delta Modulator (SAM) against a Pulse-Width Mod¬ 
ulator (PWM) shows that the SAM achieves a smaller steady-state error in a shorter 
amount of time, while using less propellant. These results, validated through numeri¬ 
cal simulations, suggest that Sigma-Delta modulation is suitable for spacecraft thruster 
control as it offers better command signal tracking. However, due to the short period 
of the SAM’s limit cycle, the propellant consumption during the steady-state increases 
at a higher rate, when compared to a PWM. To avoid this drawback, the SAM can 
be coupled with a deadzone or replaced by a lower-bandwidth modulator when the 
steady-state is reached, thus increasing the period of the limit cycle and thereby reduc¬ 
ing fuel consumption. The linearized closed-loop SAM model is found to exhibit similar 
response characteristics to both the experimental and simulated nonlinear responses, 
allowing it to be used, in conjunction with linear analysis techniques, for initial sizing 
and response estimation. 
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CHAPTER 5: 

Experimental Evaluation Methodology 


An experiment is a question which science poses to nature, and a 
measurement is the recording of nature’s answer. Before an 
experiment can be performed, it must be planned-the question to 
nature must be formulated before being posed. Before the result 
of a measurement can be used, it must be interpreted-nature’s 
answer must be understood properly. 


—Max Planck 

“The meaning and limits of exact science” [131] 


A journal version of the work presented on the Standard Test Framework is currently un¬ 
der review for publication in the Institute of Electrical and Electronics Engineers (IEEE) 
Transaction on Control Systems Technology [71], A journal version of the work pre¬ 
sented on the proposed Guidance Comparison Metric is under review for publication to 
the AIAA Journal of Guidance, Control, and Dynamics [75j. 


5.1 Introduction 

With the vast number of RT-G&C methods and emerging number of CG&C algorithms 
(a term recently coined by Lu [49j) available in literature, the challenge arises of de¬ 
termining the “best” method for a given application. While traditional RT-G&C meth¬ 
ods tend to use algebraic control laws, CG&C methods rely heavily on a model- or 
data-driven iterative process and are typically formulated using optimization-based ap¬ 
proaches [49j, [50j. The challenge then arises of selecting the appropriate G&C ap¬ 
proach that not only enables completion of the mission objectives at hand, but yields 
an acceptable level of efficiency while simultaneously meeting system-specific require¬ 
ments and constraints (e.g., algorithm reliability, computational time/responsiveness). 
Subsequently, a need arises for an equitable methodology to characterize algorithms 
in a meaningful and useful manner. 

The first half of the characterization consists of constructing relevant benchmarking 
scenarios in an effort to understand the inherent characteristics of a G&C method. 
To enable an objective evaluation (and validation) of real-time guidance and control 
algorithms, a standard framework and set of benchmarking tests must be used. The 
notion of a Standard Test Framework, along with various conceptual test scenarios 
for rendezvous and docking, to evaluate guidance and control algorithms was initially 
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proposed by Virgili-Llop et al. [52], The proposed scenarios in Section 5.2.2 build 
upon this notion and formally defines a Standard Test Framework that can be used to 
experimentally evaluate rendezvous and docking methods. 

The second part of this methodology consists of selecting and understanding the mea¬ 
surements and metrics used to characterize a given RT-G&C. Notably, a measurement 
is the quantification of some feature or attribute of the system (typically expressed as 
a single number) while a metric can be more “abstract” and is composed of one or 
more measurements [132]. When reporting on a new G&C method, the majority of 
measurements presented in literature (e.g., control effort, maneuvering time, computa¬ 
tional time) are implementation dependent, thus making it difficult to compare guidance 
methods. Therefore, it is necessary to combine these measurements in such a way 
as to form a metric that is agnostic to any specific implementation, denoted here as 
a “global" metric. To the authors’ best knowledge, no such global comparison metric 
which characterizes G&C methods has been proposed in literature. 

The remainder of this chapter will detail the proposed Standard Test Framework which 
will be utilized as the benchmarking scenarios to evaluate RT-G&C algorithms. A gen¬ 
eral formulation and interpretation of the Guidance Comparison metric is presented. To 
illustrate the extensibility of this metric, formulations for both a spacecraft and terrestrial 
vehicle are given. 


5.2 Standard Test Framework 

5.2.1 Overview 

The proposed test framework specifies the following information for each test: 

• Information available to the guidance algorithms; 

• Initial state (position, velocity, and orientation) of the Chaser and Target; 

• Docking cone corridor parameters; 

• Location, size, and movement of any obstacle(s) and associated keep-out zones 
in the workspace; and 

• Tuning objectives 

The following proposed scenarios are relevant as they evaluate the ability of the guid¬ 
ance method to handle various combinations of constraints while simultaneously pro¬ 
ducing safe trajectories and reducing propellant usage. 
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5.2.2 Test Scenarios 

A schematic of the four proposed path-constrained test cases is provided in Figure 5.1 
while the associated test parameters are listed in Table 5.1 and Table 5.2. It is worth¬ 
while to comment briefly on the selection of some the parameters. First, the location 
of both the Chaser and Target were chosen in an effort to make the best use of the 
POSEIDYN test bed’s workspace. Next, the docking cone corridor half-angle was cho¬ 
sen to be similar to the final approach cone corridor of the ISS [18]. To facilitate the 
tuning of guidance algorithms, a maximum rendezvous time is imposed only on Case 
0, which consists of an obstacle-free, straight-in approach. Note, the maximum ren¬ 
dezvous time of 180 seconds corresponds to average speed of just over 0.02 m/s - 
which is slower than the average speed of the ATV over the final approach segment 
to the ISS [18]. Additionally, prior to entering the docking corridor, the attitude of the 
Chaser must be within 10 deg^ of the final docking attitude, 9f - thus ensuring a safe 
and successful docking. Lastly, to add an additional layer of realism to each test case, 
the only information available to the guidance algorithms about each obstacle is its 
instantaneous position and velocity. 

The first test case. Case 0, is a straight-in case with only the docking cone corridor 
constraint. A time constraint of 180 seconds is imposed only on this test scenario to 
facilitate the selection of tuning parameters for an algorithm. This test case serves as 
the baseline for comparison with other cases. 

Case 1 was adopted from [52], [72] which consists of a single obstacle placed directly 
in-between the Chaser and Target. Note, this is a pathological test case for the a 
double integrator system as the Chaser will encounter the obstacle along the straight- 
line path connecting the initial and final positions. This will result result in a non-unique 
solution to guidance and control problem. Additionally, inclusion of this test case will 
allow for comparison with previously published results obtained using IDVD and Linear- 
Quadratic MFC (LQ-MPC) methods. 

Case 2 was adopted from Park et al. and incorporates two closely-spaced obstacles 
that the Chaser must navigate through (or around) in order to reach the Target [61]. 
Unlike Case 1, the obstacles are not located directly in-between the Chaser and the Tar¬ 
get, thus providing insight into the path-planning capability of the guidance and control 
algorithm. Inclusion of this test case will allow for comparison with previously published 
results obtained using both a LQ-MPC and a Nonlinear MPC (NMPC) method. 

Lastly, Case 3 is composed of an obstacle moving in a circular path that is centered 
on the straight-line connecting the initial position of the Chaser and Target. The center 
of this circular is located close to the entry of the docking cone corridor. The rotation 
rate of the obstacle was chosen such that the obstacle completes one revolution in the 

^This corresponds to the approximate half-angle of the FSS docking cone. 
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maximum amount of time allotted to complete Case 0 (180 s). This rotation rate will 
ensure the obstacle moves into and out-of the path of the Chaser multiple times over 
the duration of the maneuver. Additionally, this stressing case will provide insight into 
the ability of the algorithm to adapt to a dynamic environment. 


Case 0 Case 1 



01234 01234 

X (m) X (m) 


Case 2 Case 3 



X (m) X (m) 

Figure 5.1. Schematic of the Standard Test Framework test cases. 


5.2.3 Standard Measurements 

The proposed standard measurements associated with this framework are: control 
effort, rendezvous time, constraint handling (i.e., constraint violations), as well as the 
computational cost [52], [72]. The control effort, Ce, is defined as 

Ce= [ ||u(t)||i dt (5.1) 

Jto 
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Table 5.1. Standard Test Framework parameters. 


Parameter 

Value 

Chaser Initial Location, 

(3.5, 3.5) m 

Chaser Initial Velocity, pc 

(0, 0) m/s 

Chaser Initial Orientation, 9c 

0 deg 

Target Location, p^ 

(0.5, 0.5) m 

Target Velocity, pt 

(0, 0) m/s 

Target Orientation, 9t 

45 deg 

Target Keep-Out Radius, n koz 

rt,Koz > 0.75m 

Docking Cone Corridor Half-Angle, 13 

10 deg 

Docking Cone Corridor Length, 4 

0.75 m 

Docking Cone Corridor Entrance Attitude 

9f ± 10° 

Speed at Docking, v/ 

V/ < 0.1 m/s 

Maximum Rendezvous Time for Case 0 

180 s 

Case 3 Obstacle Center of Rotation, Pc,obs 

(2.0, 2.0) m 

Case 3 Obstacle Radius of Rotation 

0.25 m 

Case 3 Obstacle Initial Angular Position, 6*o,obs 

0 deg 

Case 3 Obstacle Angular Rotation Rate about Pc,obs 

-2 deg /s 


Table 5.2. Obstacle location and size for each test case. 


Case 

Obstacle Location (m) 

Keep-Out Radius (m) 

0 

1 

(2.5, 2.5) 

0.40 

2 

(2.0, 2.5), (3.0, 2.5) 

0.30 

3 

(2.0, 2.0) 

0.25 


where u(t) is the control input vector. The Li—norm of the control input provides a 
quantitative measure of the efficiency of the RT-G&C method and is measured in N-s. 
It is worthwhile note, other useful quantities to measure the control effort, such as 
thruster ON-time, can be derived from this measure [52], [72]. Next, the maneuvering 
time is measured as the time elapsed from enabling the Chaser onboard guidance to 
completing the maneuver (e.g., docking with the Target). Note, these two quantities 
can be used to derive the non-dimensional control effort metric, ui proposed and dis¬ 
cussed in Section 5.5.1 Lastly, the computational cost of the guidance algorithm can 
be estimated by measuring the CPU time taken to derive a control input from the inputs 
to the algorithm. 
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5.3 General Formulation of the Guidance Control Met¬ 
ric 

To fully describe a G&C algorithm, it is necessary to capture metrics spanning several 
categories including computational complexity, reliability, correctness (i.e., both numer¬ 
ical and temporal accuracy [24]), and system performance [133]. From these cate¬ 
gories, the typical measurements of interest when evaluating candidate algorithms for 
a particular application is the total (expended) control effort, computational time, and 
maneuvering time. Steinfeldt et al. proposed a framework utilizing concepts derived 
from Multidisciplinary Design in an effort to develop a parametric-based comparison 
index for various guidance and control algorithms with a common application [133]. 
This methodology, however, is currently limited to a simulation-only environment as it 
requires a Monte Carlo analysis to produce a parametric index. Additionally, bias is 
invariably introduced into the comparison index as the weighting factors for each mea¬ 
surement (indicating their relative importance) are selected through heuristics [133], 
[134]. To the author’s best knowledge, no comparison metrics have been proposed 
in literature which effectively capture the metrics of interest and compares them for 
various G&C algorithms across various system architectures. 

Drawing inspiration from the supercomputing community [134]-[136], a proposed met¬ 
ric to compare different G&C algorithms has the following qualities: 

1. Unbiased: This property ensures no one system architecture is favored over an¬ 
other architecture [134]. For instance, the metric should be invariant to different 
system architectures (e.g., different acceleration capabilities) for the same ma¬ 
neuver. As noted by Hsu et al., this property can be rather subjective [134]. 

2. Bounded: The goal of this property is to ensure the metric has either an upper or 
lower bound. This bound should also indicate either an optimal or desired value 
of the quantity. 

3. Insightful: Depending on the interpretation of the metric, its value should provide 
some insight into how well the system is performing. Additionally, an insightful 
metric can also be used to motivate design choices [134]. 

The challenge arises of selecting an appropriate metric which effectively captures the 
three primary measurements of interest, the total control effort Ce, the maneuvering time 
Atm, and the computational time of the G&C method, Tg. A scalar composite metric, 
J e M, can be formulated using either an additive form [133], 

J = \a\ce + \P\Atm + hlTg (5.2) 
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or a multiplicative form [134]-[136], 

J = (ce)“(At™)^(T,)^ (5.3) 

where a,/3,7 G M are the relative weighting factors for the measurements Ce, 
and Tg, respectively. As illustrated by these two forms, a scalar composite metric is 
clearly biased; for a given set of measurements, the value of J can change by ad¬ 
justing the relative weighting factors [133], [134]. A vector-valued metric, however, 
mitigates the introduction of bias as it implies each component of the metric is of equal 
importance [134]. Since the performance of the algorithm (e.g., expended control ef¬ 
fort and maneuvering time) are considered to be just as important as the computational 
performance of the algorithm, the use of a vector-valued metric is appropriate. 

The proposed Guidance Comparison metric has the form: 

z = [z/i, z/ 2 ]''' e (5.4) 

The first element of the Guidance Comparison metric, a time-averaged non-dimensional 
control effort, is defined as The first component of z is the non-dimensional control ef¬ 
fort defined as 

Ui = K -(5.5) 

where Cg control effort exerted by the system to complete the benchmarking maneuver 
in Atm seconds, m G M+ is the mass of the system; go is the gravitational acceleration 
constant at sea-level, 9.81 m/s^; and k g M+ is a strictly-positive factor which makes 
z/i non-dimensional. Note, the non-dimensionalizing factor, k, is a system-specific con¬ 
version factor as it is dependent upon the dimensionality of the control effort. This will 
be discussed in Section 5.5. Without loss of generality, the control effort is a system- 
specific unit of measure needed to complete a benchmark maneuver (e.g., the integral 
of the Li-norm of the control input). By having a general definition for z/i, it becomes a 
metric which can be applied across several domains. Additionally, to form an unbiased 
performance metric, any dependencies between the maneuvering time and control are 
decoupled by using them together to form a scaled time-averaged control. Scaling z/i by 
the system mass attempts to eliminate any bias due to different system designs (e.g., 
thrust-to-mass ratios), enabling a more global comparison of algorithms. This scaling 
is discussed in Section 5.5. 

The second element, iy 2 , of the comparison vector z is a non-dimensional sampling 
time 

^2 = ^ (5.6) 

-Lb 
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where Tg is the computational period of the guidance method (i.e., the guidance sam¬ 
pling period) and Th is the duration of the fastest task of the underlying real-time oper¬ 
ating system, denoted here as the Executive sampling time. The ideal value and lower 
bound for this parameter is unity. This implies the guidance algorithm is running at 
the fastest possible rate, the Executive Rate, for a given (real-time) system. Note, the 
inverse of this metric, can be viewed as a measure of the algorithm’s “real-time effi¬ 
ciency”. Additionally, this metric also provides insight into the computational complexity 
of the algorithm. A larger guidance sampling period implies the algorithm requires a 
larger amount of time to solve - thus implying greater computational complexity. 


5.4 Interpretation of the Guidance Control Metric 

Similar to benchmarking software, comparing several G&C algorithms for a given appli¬ 
cation requires the use of a standard benchmark in order to form an equitable compar¬ 
ison. For spacecraft AR&D or RPO, the benchmark of choice could be the Standard 
Test Framework proposed in Section 5.2 and [71]. To compare several guidance al¬ 
gorithms, the elements of the comparison vector, z = can be treated as a 

pair of coordinates with z/i being the abscissa and logig (z/ 2 ) being the ordinate. The 
resulting graph simultaneously captures both computational complexity and algorithm 
performance, allowing for comparison with other algorithms. It is worthwhile to note, ad¬ 
ditional insight can also be gained by plotting z/i (as the ordinate) as a function of each 
the test case for each guidance and control algorithm, as illustrated in Section 6.9. 

The ideal G&C algorithm is one that produces the optimal solution to the desired OCR 
at the fastest rate of the underlying real-time operating system (i.e., the Executive rate). 
Since both the Minimum-Fuel and Minimum-Energy optimal control problems reduce 
the amount of control effort at the expense of increased maneuvering time, the ideal 
value, ul, is as small as possible, while still completing the maneuver. Maneuvers 
which are sub-optimal will fall to the right of this boundary and take on larger values of 
z/i > z/i*, indicating a larger time-averaged control effort. Similarly, for the Minimum- 
Time optimal control problem, the maneuver time is minimized at the expense of in¬ 
creased control effort. The resulting ideal and upper-bounding value sets the maximum 
value, ul, that can be achieved for this particular objective function. Values of ui < 
indicate sub-optimal solutions which take a larger amount of time compared to the op¬ 
timal solution. Likewise, the lower bound on the non-dimensional sampling period, z/^, 
is one. Values of z /2 > 1 indicate the algorithm will be less responsive than the com¬ 
putational architecture will allow. That is, the algorithm can only respond once every 
z /2 Executive sampling periods to any perturbations and disturbances. Notional plots 
for a Minimum-Fuel/Energy and Minimum-Time optimal control problem illustrating this 
interpretation of the guidance control metric are illustrated in Figure 5.2. 
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Figure 5.2. Notional illustration on how to interpret the Guidance Compar¬ 
ison Metric. 


5.5 System-Specific Formulations 

To show the extensibility of the proposed characterization framework, this metric is 
applied to several systems. A case study focusing on the development of tradespace 
for the selection of guidance and control algorithms is presented in Chapter 7 


5.5.1 Spacecraft 

For a spacecraft, the control effort, Ce, to complete an arbitrary maneuver is defined as 
follows: 

/'(IlFWIli + lIrWIli)* (5.7) 

Jto 

A fundamental issue arises as the control effort consists of forces and torques . To 
resolve this issue, a torque can be represented by a pair of opposing thrusters of equal 
strength, F, that have a moment arm, bo, of 1 m from the center of mass of the space¬ 
craft and generate a pure rotation (i.e., a couple), as illustrated by Figure 5.3. The 
equivalent force generated by each thruster is 
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Figure 5.3. Illustration of the conversion of a torque to a force couple. 

The control effort in Equation (5.7) can now be expressed a combination of the trans¬ 
lational forces (denoted as z/i t) and an equivalent “rotational force” (denoted as 



(5.9) 


Substituting (5.9) into Equation (5.5) and taking n to be unity, the non-dimensional 
control effort, z/i, for a spacecraft becomes 



(5.10) 


where ho is 1 m. It is worthwhile to note that Equation (5.10) takes the form of a thrust- 
to-mass ratio. This is important as the 

Using the Taguchi analysis method^ [137]-[139], the effect of varying the thrust-to-mass 
and torque-to-mass ratios on vi was explored by solving for the optimal control (using 
GPOPS - II [89]) to complete Case 0 of the Standard Test Framework (described in 
Section 5.2). The Taguchi method provides a framework in which to assess the ef¬ 
fect on the output of a process by varying selected inputs using a minimal number 
of experiments [137], [138], [140]. For this study, the control variables (i.e., process 
inputs) are chosen to be the thrust-to-mass (T/m) and torque-to-inertia (r/J) ratios; 
the performance criterion (i.e., process output) was chosen to be the non-dimensional 
control effort vi defined in Equation (5.10). The magnitude of the range of the con¬ 
trol variables was selected to be 100 by surveying relevant RPO and formation flying 
spacecraft missions [7], [11], [141]-[150]. This value corresponds to the approximate 
magnitude range of thrust-to-mass ratios surveyed in Table 5.3 and illustrated in Fig¬ 
ure 5.4. In order to capture any possible nonlinearities in the performance criterion 
response, three discrete levels, { 1 / 100 , 1 , 100 }, for the control variables were used. 
Resultantly, a 19 ( 8 ^) (two factors at three levels resulting in nine experiments) Taguchi 
orthogonal array^ was used [138]. 

^See Appendix B for more information on the Taguchi analysis method. 

^See Table B.1 for the structure of the Lg(3^) orthogonal array. 
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Figure 5.4. Distribution of thrust-to-mass ratios for surveyed spacecraft. 
Adapted from [7], [11], [141]-[150]. 
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Figure 5.5. ui scaling study results showing its scaling invariance. 


To simplify the (numerical) experiments, the nominal mass and inertia were chosen to 
be 1 kg and 0.1 kg-m^, respectively; the maximum force and torque was then scaled 
according to this input. The design of experiment (DOE) was performed using both a 
minimum-propellant and minimum-energy cost functional in the optimizer. This was 
done in order to validate the results (via comparison) between the computationally 
tractable minimum-energy and (less tractable) minimum-propellant objective functions. 
The average effects of the two control variables on the performance criterion (i.e., z/i) 
was then computed as described in [137] and normalized by the average output of ui 
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across all experiments using the minimum-propellant cost function. Lastly, to assess 
the validity of the solutions, the range of values of the Hamiltonian was computed. 

As illustrated by Figure 5.5b, the range of the Hamiltonian for each cost function across 
the nine experiments is small; as such, the Hamiltonian can be considered constant 
and the resulting solution valid [76]. As illustrated by Figure 5.5a, the non-dimensional 
control effort, z/i, given in Equation (5.10), was demonstrated to be invariant to the 
scaling of the maximum rotational and translation accelerations for a relevant magni¬ 
tude range. 

5.5.2 Terrestrial Vehicles 

To illustrate the extensibility of this framework, the guidance comparison metric is for¬ 
mulated for autonomous (terrestrial) vehicles (e.g., automobiles and aircraft). A com¬ 
mon metric used to measure the control effort of an autonomous vehicle is the power 
output of the engine. Pout, [151 HI 55], 


Ce= [ Pou\{t)dt (5.11) 

Jo 


To non-dimensionalize z/i, k is chosen to be inversely proportional to the average ve¬ 
locity, V, of the vehicle throughout the maneuver. The resulting non-dimensional control 
metric is 


z/i = 


1 

V go m{tf - to) 


rtf 


Pout(t) dt dt 


(5.12) 


Scaling the work performed by the vehicle to complete the maneuver by v~^ results in 
Equation (5.12) effectively becoming equivalent to a force-to-weight ratio - resulting in 
a vi similar to that of a spacecraft, effectively turning a vehicle into a spacecraft! 


The mean power-to-weight ratio of 336 automobiles surveyed between 1975-2016 was 
found to be 0.0696 kW/kg ± 0.0143 kW/kg with a magnitude range of 2.22 [156]. The 
distribution of the surveyed automobile power-to-weight ratios are illustrated in Fig¬ 
ure 5.6. Compared to a magnitude range of the spacecraft thrust-to-mass ratios of 
113.78, the power-to-weight ratio for terrestrial vehicles is significantly smaller. Since 
the magnitude range of terrestrial vehicles is much smaller than that for spacecraft, the 
scaling will behave similar to that observed in Figure 5.5a. 


5.6 Summary 

As suggested by Max Planck, before various guidance algorithms can be compared, 
a (standard) benchmark and comparison framework must be considered [131]. The 
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proposed Standard Test Framework defines relevant scenarios that test various capa¬ 
bilities of candidate scenarios in a relevant and dynamic environment. The Standard 
Test Framework is used in Chapter 6 to evaluate and compare several guidance and 
control algorithms. In addition, the proposed global Guidance Comparison Metric pro¬ 
vides a unique capability to compare the results of different algorithms completing the 
same maneuver. This proposed metric, never seen before in literature, allows for an 


Table 5.3. Summary of surveyed RPO spacecraft design specifications. 
Adapted from [7], [11], [141]-[150]. 


Spacecraft 

Mass 

Design 

Attitude 

RPO Thruster 

Thrust-to-Mass 

(kg) 

Mission 

Control Method 

Force, (N) 

Ratio, (mm/s^) 

ASTRO 

1,089 

Formation 

Thruster 

3.5 

3.67 

(Orbital Express) 

Flying 

ATV-1 

20,252 

ISS 

Thruster 

220 

10.86 


Resupply 




Cluster 

1,184 

Science 

Spin-Stabilized 

10 

8.45 

Cygnus Orb-4 

7,492 

ISS 

Resupply 

Thruster 

31 

4.14 




3-Axis 



Hayabusa 

510 

Science 

Stabilized 

22 

43.14 




3-Axis 



Hayabusa 2 

590 

Science 

Stabilized 

22 

37.29 

HTV-2 

16,500 

ISS 

Thruster 

110 

6.67 


Resupply 




HTV-4 

16,000 

ISS 

Thruster 

120 

7.50 


Resupply 

Magnestospheric 

1,250 

Science 

Spin-Stabilized 

4.5 

3.60 

Multiscale (MMS)^ 

18 

14.40 




3-Axis 

0.8 

0.38 

OSIRIS-Rex^ 

2,110 

Science 

Stabilized 

4.5 

2.13 




22 

10.43 

Mango 

150 

Formation 

3-Axis 

1 

6.67 

(PRISMA) 

Flying 

Stabilized 

Progress^ 

7,290 

ISS 

Thruster 

9.8 

1.34 

Resupply 

Thruster 

98 

13.44 

Prox-1 

61.5 

Formation 

3-Axis 

0.05 

0.81 

Flying 

Stabilized 

Rosetta 

3,065 

Science 

Thruster 

10 

3.26 


Geometric Mean, mm/s^ 5.29 

Standard Deviation, mm/s^ 11.10 

Range,mm/s^ 42.76 

Magnitude Range 113.78 

^ Multiple thrusters used for RPO. 
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Figure 5.6. Power-to-weight ratios of surveyed automobiles. Adapted from 
[156]. 


unbiased comparison of algorithms across multiple system architectures. Due to its 
generic formulation, the Guidance Comparison Metric can be applied to virtually any 
type of mechanical system. A case study exploring one possible use of the Guidance 
Comparison Metric is illustrated in Chapter 7. 
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CHAPTER 6: 

Evaluation of Real-Time Spacecraft Guidance 
and Control Methods 


The striking difference of what from now on we will succinctly refer 
to as algorithmic control (AC) with prior practice is the absence of 
a control law ... that, before even implementation, has decided ... 
the appropriate resulting control action for each sensory input. ... 
On the contrary, AC departs from this paradigm and delegates the 
control action, including the “control law” design itself, to an online 
algorithm during execution. 

—Panagiotis Tsiotras and Mehran Mesbahi 
“Toward an algorithmic control theory” [50] 


A journal version of the experimental setup and results for APF and AAPF methods pre¬ 
sented in Section 6.3.1 is currently under review for publication in the IEEE Transaction 
on Control Systems Technology [71]. A conference version of the experimental setup 
and results for IDVD method presented in Section 6.3.2 and the LQ-MPC/RH method 
presented in Section 6.3.3 was originally presented at the 6th International Conference 
on Astrodynamics Tools and Techniques (ICATT) in Darmstad, Germany (14-17 March 
2016) [52]. A journal version of this work has been published in the CEAS Space 
Journal [45]. A conference version of the experimental setup and results for the LQ- 
MPC/DH and NMPC methods presented in Section 6.3.3 was originally presented at 
the AIAA/AAS Astrodynamics Specialist Conference 2016 in Long Beach, CA (12-16 
September 2016) [61]. Lastly, a conference version of the experimental setup and re¬ 
sults for the proposed closed-loop RRT* method was originally presented at the 27th 
AAS/AIAA Spaceflight Mechanics Conference in San Antonio, TX (February 2017) [73]. 


6.1 Introduction 

The challenge of efficiently performing a rendezvous and docking with a Target space¬ 
craft in a cluttered and dynamic environment is addressed in this chapter. This particu¬ 
lar problem incorporates the need for an algorithm to satisfy both system and path con¬ 
straints (e.g., actuator saturation limits, obstacle keep-out zone constraints) in addition 
to reducing the fuel expenditure. To accomplish this, several guidance algorithms will 
be evaluated and compared using the proposed Standard Test Framework. This frame¬ 
work, introduced in Section 5.2, consists of four test scenarios that establish a baseline 
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and evaluates the path-planning and re-planning abilities of a guidance method in both 
a static and dynamic environment. Furthermore, the framework establishes several 
measurements which can be leveraged to compare the guidance algorithms, such as 
the control effort, maneuver time, computational cost, as well as the proposed Guid¬ 
ance Comparison metric, z = [z/i, z/ 2 ]^. Utilizing a facility such as the NPS POSEIDYN 
test bed allows the algorithm to be tested in a relevant, computationally constrained 
environment in the presence of actuator response uncertainties and sensor noises. 

The remainder of this chapter is devoted to solving the rendezvous and docking prob¬ 
lem using several guidance methods which are assessed to be the state of the art in 
literature [26]. In particular, the APF, AAPF, IDVD, LQ-MPC, and a NMPC methods are 
compared to a RRT*-based method proposed in this chapter. To serve as a benchmark 
(and lower bound), the numerical optimal solution (via GPOPS — II ) or the the pro¬ 
posed test scenarios for a FSS is presented. These results, as well as those derived 
from the experimental evaluation of the above methods will be compared utilizing the 
Guidance Comparison Metric . 


6.2 Optimal Control Solutions to the Standard Test Frame¬ 
work 


6.2.1 Optimal Control Problem Formulation 

In order to establish a comprehensive benchmark as well as a lower bound for the 
non-dimension control effort z/i, the minimum-propellant optimal control problem for 
the four test scenarios associated with the Standard Test Framework was formulated 

.iT 


and solved. Defining the state vector to be x(t) = 


^x,y,9,x,y,9^ E and the 

admissible control input u(t) = [fx,fy,TY e U c the minimum-fuel optimal 
control problem for each roto-translational test case can be stated generically as: 
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Minimize: J = 

/■*/ 

Jo 

1 u(t)| 1 dt 

(6.1a) 

Subject To: x 

= 

Ax + Bu 

(6.1b) 

|u(t)| 

< 

Umax for i — 1)2,3 

(6.1c) 

x(tf) 

= 

X/ 

(6.1d) 

if 

= 

tf 

(6.1e) 

nci ■ Pdock - rid ■ r-(t) 

< 

0 

(6.1 f) 

ric2 • Pdock - ric2 ■ r{t) 

< 

0 

(6.ig) 

Ri 11^01112 

< 

0 

(6.1 h) 

where A g and B e are the state and control-input matrices, respec- 


lively, and are defined as 


A 

B 


OsxS IsxS 
OsxS OsxS 

OsxS 

diag(l/m, 1/m, 


(6.2a) 

(6.2b) 


where Osxs denotes a 3 x 3 zero matrix, Isxs denotes a 3 x 3 Identity matrix, diag(-) in¬ 
dicates a square matrix whose diagonal elements are (■), and m and J^z are the mass 
and inertial about the x axis, respectively. The desired terminal condition is denoted 
as Xf while is the maximum time allotted to complete the maneuver. The docking 
cone corridor path constraint, represented by Equations (6.1 f) and (6.1 g), is active for 
each test case whereas the obstacle keep-out zone constraint, represented by Equa¬ 
tion (6.1 h) is only active for Cases 1-3. 


To accommodate the docking cone corridor constraint, which is only active during ter¬ 
minal approach, a multi-phase OCR was formulated. The first phase consisted of the 
trajectory up to the docking cone corridor interface and includes only the obstacle keep- 
out zone (if applicable); the second phase consisted of the terminal approach inside the 
docking cone corridor. To ensure proper linkage between the two phases, denoted by 
(■)(^) for Phase 1 and for Phase 2, the following additional constraints are im¬ 
posed: 



(6.3a) 

^( 2 ) _ .( 1 ) 

''0 ~ 

(6.3b) 

\6^p-6^p\ < 10 deg 

(6.3c) 
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The constraint specified by Equation (6.3a) ensures the final state of Phase 1, is 

equal to the initial state of Phase 2, Likewise, continuity of time is enforced by 
Equation (6.3b). The third linkage constraint given by Equation (6.3c) ensures the final 
attitude of the second phase, is smaller than the final attitude of the first phase, 

0 ^. 


6.2.2 Results 

A solution for each test case for the Standard Test Framework was found via GPOPS — II 
[89]. The resulting trajectories, as well as six equally-spaced attitude snapshots, are 
illustrated in Figure 6.1. A summary of the control effort for each test scenario is tabu¬ 
lated in Table 6.1. Additionally, the points of closest approach to the obstacle keep-out 
zone(s) are annotated with a red asterisk (*) in Figure 6.1 and presented for each test 
scenario in Table 6.2. 

As illustrated by Table 6.1, the rotational control effort is two orders-of-magnitude 
smaller than the translational control effort. Likewise, the rotational component of ui 
is two orders-of-magnitude smaller and does not significantly affect the overall z/i. As a 
result, only the translational component of ui is reported since it is approximately equal 
to the overall ui. 


Table 6.1. Summary of the optimal control solutions. 



Translational 

Control Effort 
(N-s) 

Rotational 

Control Effort 
(mN-m-s) 

Translational 

^1 

(xlO-^) 

Case 0 

0.606 

7.20 

0.345 

Case 1 

0.749 

7.07 

0.426 

Case 2 

0.606 

7.20 

0.345 

Case 3 

0.616 

7.01 

0.382 


Table 6.2. Summary of closest approach for each test case. 



Closest Approach 
Distance 
(m) 

Time of 

Closest Approach 
(s) 

Case 0 

— 

— 

Case 1 

0 .000 

80.76 

Case 2 

0.058 

46.89 

Case 3 

0.000 

101.41 
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Y (m) Y (m) 



X(m) 

(a) Case 0 





X(m) 

(c) Case 2 

Figure 6.1. Trajectories of the optimal control solution for the Standard 
Test Framework test cases with six equally-spaced (in time) attitude snap¬ 
shots. 
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6.3 Formulation of Evaluated Real-Time Guidance and 
Control Methods 

6.3.1 Artificial Potential Functions (APF) 

A journal version of the guidance formulation presented in this sub-section is currently 
under review for publication in the IEEE Transaction on Control Systems Technology 
[71] as well as presented at the 2016 AAS/AIAA Spaceflight Mechanics Conference in 
Napa Valley, CA (February 2016) [32], 


Constraint Handling 

Two types of constraints are considered: a keep-out zone surrounding an obstacle and 
a docking cone corridor. These path constraints are realized via a repulsive potential, 
that creates an area of higher potential in inadmissible regions of the configuration 
space, such as exclusion-zones surround an obstacle. Two different repulsive poten¬ 
tial functions are considered: one for obstacle keep-out zones and one for the Target 
boundary constraint. A Gaussian function is used to represent a obstacle and its asso¬ 
ciated keep-out zone as is defined in the next sub-section [19], [28]. 


The second constraint type, the docking cone corridor, not only ensures a safe, collision- 
free final trajectory to the Target, but also guarantees the docking action occurs parallel 
to the docking axis [28]. The docking cone corridor is anchored at the Target dock¬ 
ing interface and defined by two parameters, the cone half-angle, {3, and the length 
of the corridor, 4, as illustrated in Figure 6.2. To enforce the docking cone constraint, 
a cardioid-like function consisting of the portions of four ellipses^ whose boundary is 
defined by the distance, d, from the center of the Target spacecraft, x^, is given in polar 
coordinates as [28], [32], 


26 iaf cos (q;) 


d{a) = < 


al cos^ (a) -I- bl sin^ (a) 
—26202 cos (a) 
al cos^ (a) -I- bl sin^ (a) 
262O3 

a/ a| cos^ (a) + 4 bl sin^ (a) 
26103 

[ a/ a| cos^ (a) + 462 sin^ (a) ’ 


a 


e [0, f) 

« e [f, tt) 

« e [tt, f) 
a E [^, 2%) 


(6.4) 


‘'When bl = 62, the resulting cardiod is composed of three distinct ellipses 
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where the parameters 01 , 02 , 03 , 61 , 62 , the distance d, and the angle a are defined in 
Figure 6.2. To satisfy docking cone corridor constraint, the semi-major axes oi and 02 
are chosen such that the Target keep-out zone constraint intersects with the docking 
cone corridor at d{7i/2 ±(3) = 4 / cos {(3). Therefore, the length of the semi-major axes, 
oi and 02 which satisfy the docking cone corridor constraint: 



(6.5a) 



(6.5b) 


where don is the distance from the Target center of mass to the docking cone interface, 
«! = 7r/2 - 13 and 02 = t^I‘^ + 13- It is worthwhile to note, this approach can also be 
used to realize an asymmetric keep-out zone around the Target spacecraft due to both 
static and dynamic appendages such as antennae, articulated solar arrays, and robotic 
manipulators. This is done by selecting values for 61 , 62 , and 03 which encompass the 
appendage(s) and their respective workspaces [32]. 



Figure 6.2. Notional Target keep-out zone and docking cone path con¬ 
straints (blue) with overlaid docking cone constraint (red) in the Target 
body frame. 


Artificial Potential Function Formulation 

The Artificial Potential Function (APF) method utilizes the gradient of a potential field, 
which is composed of both attractive and repulsive potentials, to derive the necessary 
control inputs to reach the desired or goal position. The following development of the 
APF guidance method is adapted from [19], [28], [32]. 
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The attractive potential, (t)a, establishes a global minimum in the workspace at the de¬ 
sired (terminal) configuration, pj = [x/, y/, 6*/]^. The attractive potential was chosen to 
be a quadratic function given as 


<t>a = Y^lfQarcf ( 6 . 6 ) 

where r^f = Pc-P/ is the relative difference in the Chaser configuration, p^ = [x, y, 9]^, 
and the desired configuration, ka e M+ is a strictly positive constant, and G 
is a symmetric, positive-definite, attractive potential shaping matrix. In order to avoid 
certain portions of the configuration space, the repulsive potential, 0^, creates an area 
of higher potential effectively “pushing” the Chaser away from these regions [19]. Two 
different repulsive potential functions are considered — one for an obstacle keep-out 
zone and one for the Target boundary constraint. Keep-out zones centered around an 
obstacle are represented using a Gaussian function [19], [20], [28], 


No 

= exp 

i=l 




(Ti 


(6.7) 


where No is the number of obstacles, Vca is the relative position of the Chaser with 
respect to the center point of obstacle, o-j and are the width and height parameters 
for the obstacle, and Nj g is a positive definite shaping matrix for the 
repulsive potential. For this application, the repulsive potential shaping matrix is taken 
to be the Identity matrix, Isxs, implying a circular obstacle keep-out zone. 


In order to ensure the Chaser appropriately enters the docking cone constraint and re¬ 
mains clear of any appendages attached to the Target, the following repulsive potential 
function was selected as [32], 


kr rl/Q&rc/ 

2 exp - 1 ] 


( 6 . 8 ) 


where kr G M+, Q,b,Pb G are positive definite shaping matrices for the Target 
boundary constraint, and is the relative position of the Chaser with respect to the 
Target boundary constraint inertial position, x?,. Note, explicit knowledge of x;, is not 
necessary as the relative position vector Vcb can be expressed as 


Tob = d{a) 



- ret 


(6.9) 


where is the relative position between the Chaser and Target, d{a) is given in Equa¬ 
tion (6.4) and x* is the Target state. From a geometric prospective, the shaping matrices 
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Q and P can be viewed as shaping the height and width of the potential function, re¬ 
spectively. The resulting total potential function, (ptot, is defined as the superposition of 
the attractive and all repulsive potential functions: 


No 

0tOt = (f>a + 4>b + 4>ri 
i=l 


( 6 . 10 ) 


Since finer control is desired when in close proximity of the Target and desired terminal 
state, a continuous feedback control law was chosen and is given as [32], [35], 

u(Pc,P/,Pc) = -B^^Ka (pc-k Va:0tOt) (6.11) 

where B = diag[Vm, Vm, ^/h], e is a positive gain matrix and V^^tot is 
given as 


Va:0tot = haQa^cf + K exp [l - r]^PbVcb] {Qb^^cf “ (rJ^Q^rc/) Pbrcb) 


No 


-E 


2^|Ji 


2=1 


(Ti 


exp 


(Ji 


Nr 

^ 2 ^ CO 


( 6 . 12 ) 


Stability analysis of the APF control law showing Xc approaches x/ asymptotically can 
be shown via a Lyapunov analysis and is presented in [28], [35]. 


Adaptive Artificial Potential Function Formulation 

The Adaptive APF (AAPF) method [19], [20] was developed in an effort to reduce fuel 
consumption while still providing effective and robust real-time obstacle avoidance. To 
achieve this, a time-varying attractive potential shaping matrix, Qa(t), is updated ac¬ 
cording to an adaptive update law in order to follow a given reference trajectory. The 
upper bound on the optimality the AAPF method can achieve is bounded by the op¬ 
timality of the reference trajectory. Note, the following development of the adaptive 
update laws are adapted from [19], [20]. 

To begin, the reference trajectory was generated by considering a rest-to-rest, obstacle- 
free, fixed-time, minimum-fuel OCP for a double-integrator system with a bounded con¬ 
trol input. The resulting optimal control can be described as a bang-coast-bang maneu¬ 
ver, whose desired control input direction is anti-parallel to the r-c/ [76]. Therefore, the 
resulting reference trajectory follows the path described by Pref(t) = -refit) and a ref¬ 
erence velocity, Pref(t) = -VrefTefit). Notably, for this application, the reference speed, 
Wref, can be chosen arbitrarily within the capabilities of the system. For this study, the 
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reference speed was selected to be constant by specifying a final time as described in 
more detail in Section 6.6.1. 


Since the translational and rotational motions of a spacecraft can be considered uncou¬ 
pled, two independent adaptive laws can be developed. Consider first the translational 
motion of the Chaser with the configuration variable being p = [x.yY. The Cholesky 
factorization Qa(t) = R^(t)R(t), where R(t) is the upper triangular matrix [19], [20], 


R(t) 


Pllit) pl2{t) 

0 p22{t) 


(6.13) 


is utilized to enforce the symmetric positive-definite condition of the time-varying attrac¬ 
tive potential shaping matrix. 

The weighting factors composing R(t) are then chosen such that increased efficiency 
is achieved while using the same continuous control law given in Equation (6.11) [19], 
[20]. This is achieved by first defining an error function as the difference between the 
reference velocity, pref and the negative gradient of the attractive potential 


6 Pref ( ^ x4^a} 

= -VreiTT^^ + kaBjRr. 


(6.14) 


c/ 


ikc/lla 

The time derivative of the error function for a stationary terminal configuration is 


e = -Vref 


Pc r^f^cfPc 


Ikc/ll: 


-I- ka ( R'''R -I- R'''R ) p 


(6.15) 


Letting p = [pn, pi 2 , p 22 ]^ the second term of Equation (6.15) can be parameterized 
as. 


ka R'''R -I- R'''R ) Xc = /CaKp 


where K g 


d2x3 


IS 


K = 

and whose elements are given as 


kii ki2 ki3 
^21 ^22 ^23 


ku = 2piix -f P 12 P, ki 2 = Piip, /cis = 0 

k 21 = P 12 X, k22 = PllX + 2pi2P, k23 = 2p22y 


where x = Xc — Xf and y = Uc — yj- 
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In order for the error to converge to zero asymptotically, it is desired to have e = -ke 
where /c G M+ [19], [20]. The resultant adaptive update law for specifying the attractive 
potential weights is given as follows: 

p = t;,efK (KK^)-' (- K (KK^)-' {kaK^Rp, - ke) (6.16) 

Fc/lls J 

which will drive e —)■ 0 asymptotically. 


Following a similar procedure, an adaptive update law for the rotational motion of the 
Chaser can be derived. Since the Chaser rotational motion is constrained to a single 
axis of rotation, R = The error function is given as 


60 — CJref — ( — 

9r - 9f 


— 


\0c-0^ 


f\ 


+ kaphi^c - Of) 


(6.17) 


where Wref is the reference angular speed. The resulting adaptive update law that will 
drive -)■ 0 asymptotically is given as follows: 


P33 — i^ref 


{\0c - 0f\-^ - {9, - 9f)^\\9, - 9f\-^) kapl,9, - kee 


‘^kaP33{9c — 9f) 


‘^kaP33{9c — 9f) 


(6.18) 


The adaptive update law given in Equation (6.16) and Equation (6.18) has several 
critical implications [20]. First, a nonzero initial value for p is required for an initial 
positive-definite shaping matrix. Secondly, the matrix K is non-singular for all x ^ xj. 
As the Chaser approaches its terminal state, the matrix K becomes ill-conditioned 
causing the adaptive estimates to diverge and tend to infinity [20]. Resultantly, the 
adaptive estimates are bounded by first defining a convex set for each estimate [20] 

Ai = {peR \ p~ < p<pi} 

given an upper and lower bound p'^,p~ for each estimate and letting 

P = [Pll-, Pl2, Pl3, P22, P23, PSS]"'' 
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be the augmented vector of adaptive estimates. Next, a projection algorithm is used to 
ensure the adaptive estimates remain bounded [20], 


P = Proj (p) 


Pi if Pi ^ ^i 

Pi if Pi = Pi and p* > 0 

Pi if pi = pt and p* < 0 

0 otherwise 


(6.19) 


Bounding the adaptive estimates allows the algorithm to be used for close-proximity 
and rendezvous maneuvers. 

6.3.2 Inverse Dynamics in the Virtual Domain (IDVD) 

A journal version of the guidance formulation presented in this sub-section has been 
accepted for publication in the CEAS Space Journal [72] as well as presented at the 
6 th International Conference on Astrodynamics Tools and Techniques (ICATT) in Darm- 
stad, Germany (14-17 March 2016) [52]. 

The IDVD method, a quasi-optimal, direct guidance method, solves the optimal control 
by parameterizing both the trajectory and time using basis functions [54], [72], [157]. 
The resulting parameter optimization determines the set number of coefficients while 
minimizing a user-specified cost functional and satisfying system and path constraints. 
The following development of the IDVD method is adapted from [52], [72]. 

The polynomials describing the time are defined as a function of a virtual time, Kt e 
[0,Ktj] as [52], [72], 


tx{n) = (6.20a) 

i=0 

nty 

= (6.20b) 

i=0 

where da^ and dbi are the coefficients for the polynomial describing time in the respec¬ 
tively X and y directions; nt^ and rity are the number of coefficients for the polynomials 
describing the 4 and ty trajectories, respectively. Note, no constraints are placed on nt^ 
and nty being the same. Likewise, the polynomials describing the x and y trajectories 
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are defined as [52], [72], 


'^x 

x{k) = ( 6 . 21 a ) 

i=0 

riy 

yif^) = (6.21b) 

i=0 

where a* and 6* are the coefficients for the x and y trajectories, respectively; and Uy 
are the number of coefficients for the polynomials describing the x and y trajectories re¬ 
spectively. Similar to the time polynomials, and Uy need not be the same. Thus, the 
OCP is converted into a parameter optimization where the coefficients a^, 6*, da^, and 
must be determined. 

Before continuing, it is important to note the three constraints that are placed on the 
time polynomials [52], [72]. The first constraint is the placed on the terminal condition 
of the time polynomials, specifically, the polynomial final time, is less than the 

specified maximum time tf. 

< tf ( 6 . 22 ) 

The second constraint only applies when different time polynomials are used. In this 
case, dh^ is chosen such that txintj) = ty{K,tj). The last constraint is placed on the 
derivative of the time polynomials such that it must be strictly positive over the entire 
range of the virtual times 


= '^idaiKd~^ > 0 V K 

i=0 

nty 

t'y{Kt) = > 0 V K 

i=0 


(6.23a) 

(6.23b) 
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Next, the polynomials describing the time derivatives of the x and y trajectories with 
respect to the time polynomials are given as follows [52], [72]: 


X 




i=0 


X 


X 


X = 


n 

X = 


it'.? it'.? ^ 

- l)aiK\ 


i=0 




i-2 

t 


i=0 


. y 




i=0 


.. ^ ji _ y_.>i 

^ ifx {t'r ' 




// 

y = 


- l)biKl 2 


i=0 

nty 

y''i^? = ?dhA~‘^ 

i=0 


(6.24a) 

(6.24b) 

(6.24c) 

(6.24d) 

(6.24e) 


It is worthwhile to note, the number of polynomial coefficients can be reduced given the 
(estimated) state of the Chaser and the desired terminal state [52], [72] - resulting in a 
smaller parameter optimization process. Specifically, the initial and terminus conditions 
of the trajectory-describing polynomials for the rendezvous and docking scenario are 

x(0) = Xc, i(0) = ic ff. PCX 

Atf) = X/, ±itf) = 0, ±{tf) = 0 

As a result of these initial and terminal constraints, the minimum order for the trajectory 
polynomials while still meeting system and path constraints is n.^y > 5. Additionally, 
the minimum order of the time polynomial is > 1- 

Selecting the cost functional to be the Li-norm of the control inputs, the resulting 
instantiation of the IDVD method is a minimum-propellant, bounded-time optimal control 
problem: 
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/■*/ 

Minimize: J = 

Jo 

||u 

Will dt 

(6.26a) 

Subject To: x 

= 

Ax -L Bu 

(6.26b) 

\u{t)i\ 

< 

Umax 

(6.26c) 

\u{t)2\ 

< 

Umax 

(6.26d) 

^(tf) 


= ^f 

(6.26e) 


< 

tf 

(6.26f) 

Jai,bi 

> 

0 

(6.26g) 


= 


(6.26h) 

Ri ll^ci||2 

< 

0 

(6.26i) 

nci ■ Pdock - rid ■ r-(t) 

< 

0 

(6.26j) 

nc2 • Pdock - ric2 ■ r(t) 

< 

0 

(6.26k) 


where x = [x,?/]^; u = [fx,fyY', tf is the maximum time allowed to complete the 
maneuver; and r^o is the relative position of the Chaser with respect to the obstacle. The 
docking cone corridor constraint is implemented as two hyperplanes defined by a point 
P(.) with normal vectors n(.) while the obstacle keep-out zone is kept as a nonlinear, 
non-convex path constraint. To solve the resulting nonlinear programming problem, the 
open-source Interior Point OPTimizer (IPOPT) solver [158] is implemented onboard the 
FSS at a rate of 5 Hz in order to produce a feedback action. It is important to note, the 
constraints specified in Equation (6.26) are only satisfied at the equally-spaced discrete 
nodes considered by the IPOPT solver. 

6.3.3 Model Predictive Control (MPC) 

A journal version of the guidance formulation presented in this sub-section has been 
published in the CEAS Space Journal [72] as well as presented at the 6th International 
Conference on Astrodynamics Tools and Techniques (ICATT) in Darmstad, Germany 
(14-17 March 2016) [52] and the 2016 AIAA/AAS Astrodynamics Specialist Conference 
in Long Beach, CA [61]. 

While based primarily on the LQR, the MPC framework is solves the constrained 
optimal control problem over a finite horizon subject to both system and path con¬ 
straints [55], [72] - compared to the IDVD method which optimizes over the entire 
trajectory . Feedback action is obtained by implementing the first control input and 
resolving the problem using the current estimated state of the system [52], [55], [61], 
[72]. The following development of the LQ-MPC and NMPC formulations are adapted 
from [52], [61], [72]. The MPC guidance method, as implemented, is based on the 


101 



following optimal control problem: 


Minimize: J = (xat — x/)^ P (xat — x/) 

N-l 



-L ^ (xfc+i -XfY Q (xfc+j 

^ _n 

-X/) 



t - KJ 

N-l 





i=0 


(6.27a) 

Subject To: 

^fc+1 

Axfc + Bufc 


(6.27b) 


VI 

Umax fOI" i — 1)2,. 


(6.27c) 

Ri - 1 

1 r ■ II ^ 

1 ^ C2 112 — 

0 


(6.27d) 

■ Pdock rid 

VI 

0 


(6.27e) 

ric2 ■ Pdock fic2 

■r(t) < 

0 


(6.27f) 


where state and control vectors at sampling instant k are denoted by x^ g and 
Ufc G respectively; A g and B g are the respective discrete state 

and control input matrices; r = [x,yY is the position of the Chaser; x/ is the terminal 
state; and N is the length of the horizon. The matrices Q G and R G 

define the positive definite relative weighting on the path and control effort, respectively, 
while P G is the solution to the discrete Riccati equation [55]. 

Two instantiations of the MPC method are considered, a Linear-Quadratic MPC (LQ-MPC) 
and a Nonlinear MPC (NMPC). Compared to the nonlinear variant, the LQ-based ap¬ 
proach requires both linear dynamics as well as linear inequality constraints [72]. The 
nonlinear obstacle keep-out zone constraint is linearized through two methods, a rotat¬ 
ing hyperplane method [52], [61], [72] as well as a novel dual-hyperplane method [61]. 
Compared to the LQ-based method, NMPC directly handle both nonlinear dynam¬ 
ics and constraints. Additionally, both implementations of the MPC method handle 
the docking cone corridor constraint in a similar fashion to the IDVD method. As 
implemented, the LQ-MPC method uses a publicly-available quadratic-programming 
solver [159] and resolves the optimal control problem every 3-5 s while NMPC uses 
IPQPT [158] and resolves the optimal control problem every 3 s. 


6.4 Closed-Loop Rapidly-Expanding Random Trees 

A conference version of this work was presented at the 2017 Spaceflight Mechanics 
Conference in San Antonio, TX (February 2017) [73]. 
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6.4.1 Overview of Sampling-Based Methods 

The RRT Algorithm 

The Rapidly-Exploring Random Tree, or RRT, algorithm was introduced by LaValle in 
1998 [160] as a kinodynamic path-planning method to quickly search high-dimensional, 
non-convex spaces for a feasible trajectory [64], [66]. This sampling-based path plan¬ 
ner addressed the need for a general-purpose path-planner that scales well with the 
dimension of the problem. 

Before continuing, it is worthwhile to note some key concepts and nomenclature. For 
the planar^ translational guidance problem, the configuration space is denoted as X c 
and the configuration variable be denoted as q g A tree, denoted as T, has 
a single root node go that is its initial configuration [161] and contains no cycles [161]. 
That is, each node, not including the root, has only one incoming edge, termed the 
parent node and denoted as gpar, but can be the parent to numerous other nodes, 
termed child nodes. This results in a unique path from any node belonging to the tree, 
denoted as g* G T, to the root node. Lastly, a newly added node, denoted as gnew is 
commonly referred to as a child node. [161]. 

The RRT algorithm, illustrated in Algorithms 1 and 2, begins by initializing a tree, T 
with root node, go, which is the initial configuration of the system, in this application 
the Chaser, at the start of the algorithm. The objective of the RRT algorithm is find 
a feasible, if one exists, from go to the goal configuration, g/ g X given any fixed 
obstacles, Xobs e X [64], [66], [160]. 

The tree is first initialized using the Initialize_Tree(-) function where the algorithm 
attempts to connect the root node to the goal node with a collision free path using a 
user-defined, system-specific steering function Steer(-) [160]. If such a connection is 
possible, the goal node is added to the tree and the algorithm is completed. In the 
event this is not possible, the obstacle-free configuration space, 

Xfree = X \ Xobs (6.28) 

is uniformly sampled via the Sample(-) function in line 3 to produce a random config¬ 
uration, grand- Given some user-selected metric (e.g.. Euclidean distance), the nearest 
node, gnearest> to grand 'S determined using the Nearest(-) function. For the RRT algo¬ 
rithm, the parent node, gpar, of grand is chosen to be gnearest [160]. The control input, u*, 
required to take the system from gnearest to grand is determined via Steer (gnearest, grand) 
function (line 3). If a collision-free path between gnearest and gnew exists, the resulting 

'' Note, for the three-dimensional case, the configuration space is denoted as X c config¬ 

uration variable be denoted as q e 
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configuration from applying the steering law to the system, gnew, is added to the tree; 
otherwise, the tree reamins unchanged and grand is discarded. 

This process is completed until the maximum number of iterations is obtained [160]. 


Algorithm 1 The RRT Algorithm [160]. 

1: T ^ Initialize_Tree(go) 


2 : for * = 1 doiV,^,„ 


3: grand ^ Sample(X, Xobs) 


4: T ^ Extend(T, grand) 

> gnearest IS the parent node of gnew 

5: end for 



Algorithm 2 The RRT Extend function [160]. 

1 

function EXTEND_RRT(T, grand) 


2 

ijnearest ^ Nearest(7I, Xrand) 


3 

(O'new) rij) i Stee2r(gjiearest) Q'rand) 

> gnearest IS the parent node of gnew 

4 

if 0bstacle_Free(gnearest,gnew) then 


5 

AddNode(7~, gnewr ^’nearest) ^i) 


6 

end if 


7 

return T 

> Return updated tree 

8 

end function 



It is worthwhile to note, the fidelity of collision-checking is problem specific. Collision¬ 
checking can range from a simple point-mass model to three-dimensional collision de¬ 
tection algorithms [162]-[165]. 

This method has several key advantages including excellent exploration of the search 
space, probabilistic completeness, and computational simplicity. Additionally, due to the 
construction methodology of the tree, the RRT algorithm tends to, on average, partition 
largely unexplored areas into smaller regions [64], [160]. However, the resulting path 
from the root to the goal the algorithm converges to will never be optimal [65]. In terms 
of computational complexity, the RRT algorithm is bounded by the computational effi¬ 
ciency of the nearest neighbor search [64], [66], [160]. Karaman and Frazzoli showed 
the computational complexity of the RRT algorithm at each iteration is at least on the 
order C>(log (iV^J, where Ny^ is the number of nodes at the end of the i*'* iteration [65]. 

This framework was later extended by Frazzoli et al. by exploring a new technique to 
connect the newly sampled nodes to the tree and on a method to handle moving ob¬ 
stacles. Rather than connecting the sampled node to its nearest neighbor, Frazzoli 
proposed choosing the parent based on an optimization heuristic [161]. Given an opti¬ 
mal steering law, each node in the tree, g*, is sorted based on the obstacle-free cost to 
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from qi to gnew plus the cost accumulated at node qi. Starting with the lowest obstacle- 
free cost node, the trajectory between gnew and each node in the tree is checked to 
ensure the trajectory is collision-free for some duration of time [161]. The parent node 
is selected as the first node which satisfies this condition. This method was then ap¬ 
plied to several systems including the path planning for a small autonomous helicopter 
- which was the motivation for the development of this algorithm - in addition to the 
attitude slew planning for a spacecraft with multiple constraints [83], [161]. 


The RRT* Algorithm 

Until Karaman and Frazzoli introduced an optimality extension to the RRT framework, 
no RRT algorithm was shown in literature to converge to an optimal solution. This 
extension initially introduced in 2010, called RRT* [65], [166], fuses the asymptotic 
optimality of rapidly-exploring random graphs with the tree-like structure of the RRT 
algorithm. The RRT* algorithm achieves the improved performance through the way it 
extends the tree given a new node. Specifically, the parent node is chosen from a set 
of fc-nearest nodes, denoted as Xnear, that minimizes a user-defined cost functional, 
as illustrated in Algorithm 3. 

It is worthwhile to note, there are typically three different costs are associated with each 
node {xi) [44]. First, the total cost is the cost accumulated from the root go to node g* 
and is represented by the function Total_Cost(gj). The incremental cost is the cost 
associated with going from the parent of node g* to node g*, represented by the function 
Delta_Cost(gj). Lastly, the cost-to-go is a lower-bound on the cost associated with 
moving the system from the current node to goal configuration, g/ and is typically taken 
to be the (unconstrained) Euclidean distance [44]. 

In an attempt to further optimize the tree by reducing the overall cost of the result¬ 
ing path, the cost of going from the newly added node gnew to each nearby node 
gnear ^ Xnear is evaluated. When the total cost from the root node go to the node gnear 
(Total_Cost(giiear)) is found to be greater than the total cost from go to the gnew plus the 
incremental cost go to gnear (Delta_Cost(gnear)), gnew is chosen to be the new parent 
node of gnear- When this occurs, the tree is said to be rewired [44], [65], [166], [167]. 
The rewiring process is illustrated in Algorithm 5. 

Example path planning solutions of the RRT* algorithm are presented in Figure 6.3 
utilizing a minimum fuel (Li—norm) cost function, 150 nodes, for various numbers of 
nearest neighbors, k. Note, for k = I, the RRT* algorithm becomes the RRT algorithm 
as it does not consider any additional nodes other than the closest. It can be observed, 
as the k increases, the resulting graph becomes more directed from the start point 
towards the goal and more efficient paths are found. Increasing the value of k from one 
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Algorithm 3 The RRT* Extend function. 

1: function EXTEND_RRT*(T, Xrand) 

2: gnearest ^ Nearest (T, grand) 

2- (gnew) Steer(^jjgarest) grand) 

4: if Obstacle_Free(gnearest, gnaw) then 

5: Xnear ^ Near(T, gnaw, k) 

6: (gmin, Uj) 4- ChOOSe_Parent(Xnear, gnaarest, gnaw, Ui) 

7: T AddNode(T, gmin, gnaw, ^i) 

8. i Rewire(7~, Xjjear, gmin, gnaw) 

9: end if 

10: return T 

11: end function 


Algorithm 4 The RRT* Choose.Parent function. 

1: function CHQQSE_PARENT(Xnear, gnearest, gnew) 

2: gmin gnearest ^ InitieNze gmin 3nd Cmin 

3: Cmin ^ Total_Cost(a;nearast) + Delta.Cost(x^ew) 

4: for all gnear ^ Xnear do > Select parent which minimizes the cost 

5: (g', u') ^ Steer(giiaar, gnaw) 

6: if Obstacle_Free(a;', Xnew) then 

7: c' ^ Total_Cost(gnear) + Delta_Cost(g') 

8 : if d < Cmin then 

9: Cmin ^ d ; gmin ^ g' i Umin ^ u' > Update best parent node 

10: end if 

11: end if 

12 : end for 

13: return gmin, Umin 

14: end function 


to five results in 77.1% difference between the RRT and RRT* methods; letting k be 25 
results in a 142.6% difference between the RRT and RRT* algorithms. 


Extensions and Applications of the RRT* Method 

In an effort to increase the convergence rate, Qureshi et al. modified the RRT* algorithm 
by biasing the randomly sampled points with the gradient of an APFs by placing an 
attractive potential at the goal point and repulsive potentials on the obstacles [67]. This 
method was shown to reduce the time, number of iterations, and cost associated with 
the first solution that was found. A similar approach proposed by Pharpatara et al. 
creates a convex cone from which the randomly sampled points are drawn [68]. This 
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Algorithm 5 The RRT* Rewire function. 

1 : function REVIRE(mathcalT, Xnear, a;min, a;new) 

2: for Xnear ^ Xpear \ {^^min} dO 

3: (x', u') Steer(a;new, a;near) 

4: if Obstacle_Free(a;new, then 

5: if Total_Cost(a;new) + Delta_Cost(a;') < Total_Cost(a;near) then 

6 : T ^ ReCOnneCt(T, Xnew, a^near) 

7: end if 

8 : end if 

9: end for 

10: return T 

1 1: end function 


method was found to reduce the average number of iterations by half (compared to the 
RRT* method); the average path length was reduced by approximately 15.3%. Several 
other extensions to the RRT* algorithm are detailed in [168]. 

The RRT* algorithm has been applied to several real-time systems. Karaman et al. 
extended this algorithm for use with autonomous wheeled robot by first selecting and 
storing part of the generated path, called the committed trajectory, at the end of every 
planning cycle [44]. While the robot is executing the committed trajectory, the algorithm 
re-roots the tree at the end of the committed trajectory and continues to improve the 
uncommitted portion of the path. When the committed trajectory is completed, the 
algorithm selects a new committed trajectory and repeats until the maneuver is finished. 
Additionally, this algorithm was applied to the Talos vehicle, MIT’s entry to the DARPA 
Urban Challenge [42], [169j. In this application, the connection between existing nodes 
and a newly-sampled node is performed through a nonlinear, closed-loop simulation. 
This closed-loop simulation ensures vehicle stability throughout the challenge. 

it is worthwhile to note, the systems discussed thus far utilized the RRT-based method 
strictly as a path planner. In each implementation discussed, there was a lower-level 
controller which would guides the system from one configuration to the next. It is the 
intention of this work to utilize a RRT*-based method coupled with Harmonic Potential 
Functions as a feedback controller for the translational guidance and control of a space¬ 
craft in the presence of (a) dynamic obstacle(s). The difficultly arises as a spacecraft is 
not a kinematic vehicle (i.e., not controlled by its velocity) and must manage and take 
into account its linear momentum to ensure collision avoidance. 
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(a) fc = 1, Cost:0.746 m/s 




(c) k = 10, Cost:0.150 m/s 



(d) k = 25, Cost:0.125 m/s 


Figure 6.3. Example path planning solutions using the RRT* algorithm for 
various numbers of nearest neighbors, k, and 150 nodes. 


6.4.2 Harmonic Potential Functions (HPF) Applied to the RRT 

In this work, Harmonic Potential Functions (HPFs) are utilized to generate a Naviga¬ 
tion Function which biases the generated random states towards the goal point and 
obstacle-free space. Paths are generated by following the negative gradient of the 
Navigation Function - which is commonly referred to as the velocity potential [36], [39], 
[40]. To achieve this, elementary solutions to the Laplace equation, = 0, are con¬ 
sidered [36], [38], [39]. While the primary development is for two dimensions, exten¬ 
sion to three dimensions is trivial through the use of spherical coordinates [36]. Thus, 
only elementary elements which can be extended to three-dimensions are considered. 
Since these elementary solutions have strong ties to incompressible fluid flow\ the 

^ It is worthwhile to note, solutions to Laplace’s equation shares strong ties to Electrostatics [170]. 
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subsequent development of the potential field leverages nomenclature and notation 
from Aerodynamics adopted from Reference [36]. 

To formulate the underlying navigation function (at each sampling instance), elements 
are added to the workspace incrementally, exploiting the linearity of the Laplace equa¬ 
tion and its solutions [36], thus simplifying the selection of the strengths associated with 
each element. 

First, a uniform flow that is directed towards the Chaser spacecraft from the goal posi¬ 
tion is added to the workspace [39], [40], [43], 

0uf = 'yref ||rc/|l 2 Cos(a/) 

/ f \ \ ' f w (6.29) 

= Uref (a; COS [af) + ysm [af)) 

where Wref is the reference speed of the chaser; r^/ is the relative position Chaser with 
respect to the desired terminal position; and «/ is the angle between the Chaser and 
the terminal position given as 


af = tan ^ 


f y-yf \ 

\X-Xf) 


The resulting velocity potential due to the uniform flow is 


^ 0a;,uf — '^ref 


COS (a/) 

sin {af) 


(6.30) 


The resulting potential field, illustrated in Figure 6.4, has zero potential at the goal 
position and increases radially-outward form the goal position. From Equation (6.30), 
the uniform flow will bring the Chaser from any point in the potential field to the desired 
terminal position — but it does not take into account any path constraints such as 
the docking cone corridor constraint or obstacle keep-out zone. The inclusion of the 
uniform flow eliminates numerical noise due to the small gradients associated with 
source, sink, and doublet elements at large distances [39], [40]. As it will be illustrated 
shortly for the case of a doublet, the radial component of the velocity potential of these 
elements in a polar coordinate system is inversely proportional to the distance r to 
the [36], 

V0r = ^ OC (6.31) 

dr ||r||2 


Next, to account for the docking cone corridor path constraint, a doublet - the result¬ 
ing elementary solution arising from a source and sink of equal strength which are 
made to approach each other [36] - is added to the workspace (that already contains 
the uniform) its central influx (anti-efflux direction, see Figure 6.5) pointed towards the 
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(a) 3-D Representation of the Uniform Flow Potential (b) 2-D Uniform Flow Potential Contours and Gradient 


Figure 6.4. Elementary 2-D uniform flow centered at the goal position. 


goal position and parallel to the docking axis. This creates a global minimum in the 
workspace in the vicinity of the goal position at -cx). Note, as long as the goal position 
is sufficiently close to the Target position, this is a sufficient approximation of locating 
the global minimum at the goal position as the HPF is not being directly used for nav¬ 
igation, but rather to influence the random selection of a configuration. Additionally, 
this approximation allows the use of a doublet to enforce the docking cone corridor 
path constraint in a similar fashion the cardioid function enfroces this constraint for the 
APF-based guidance methods (illustrated in Figure 6.2) [32], [71]. 

The general expression for the potential due to a doublet is 

<Pd = Tj—^ COS (7rfJ (6.32) 

IPcdi II2 

where Vdi is the relative position of the Chaser in the potential field with respect to the 
doublet; Kd^ e M is the strength of the doublet; 7 ^- is the angle between the Chaser 
and the i*'* doublet given as 


= tan 


f y-VdA 

\x-XdJ 


Note, a positive value of Kd indicates the flow proceeds outward from source-side of the 
doublet in the negative direction [36]. For simplicity, the flow direction into (influx) or out 
of (efflux) is with respect to the negative gradient of the doublet potential, as this is the 
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trajectory the Chaser will follow. That is, the directions of each flux follows the direction 
of the negative gradient. The resulting potential field and velocity potential field of a 
doublet is illustrated in Figure 6.5. 

The target doublet strength, ^or, is selected by first considering the total potential of the 
field given the two elements added thus far: 


0 = t'ref l|rc/|l2Cos(ag) + 


^tar 

Ikc/lla 


cos ( 7 t) 


(6.33) 


where Vcf is the relative position of the Chaser with respect to the desired terminal 
position. Applying the von Neumann boundary condition - which specifies the flow 
must be tangential to the boundary of a solid obstacle [40] - to the target doublet 
implies the radial velocity potential along the boundary of the target doublet is zero, 

Vr = d(l)/dr = 0, yields: 


Vr = Vref cos (Og) 


^tar 



cos ( 7 t) = 0 


(6.34) 


Note, at the boundary of the target doublet, ||r ||2 = R, where R is the radius of the 
equipotential needed to realize the docking cone corridor. Treating each equipotential 
lobe of the doublet as an ellipse (with similar semi-major and semi-minor axes), the 
necessary radius to satisfy the docking cone corridor can be found by applying Equa¬ 
tions (6.5a) and (6.5b): 


R 


+ r^off 

2 cos (/3) cos {6t) 


(6.35) 


given the docking cone half-angle, (3, corridor length,4, and offset doff and the attitude 
of the docking axis, 4. Note, for the potential field consisting of these only the uniform 
flow centered on the terminal position and the doublet centered on the Target, ag = -'ft- 
Substituting Equation (6.35) into Equation (6.34) and solving for target doublet strength 
yields: 


_ f 4 + dow \ 
\2 cos (/3) cos (4) / 


(6.36) 


Lastly, doublets can be added throughout the workspace to represent circular obstacle 
keep-out zones with radius Ri [38], [40], [43]. To enable greater safety, each doublet is 
free to rotate such that the negative velocity potential (i.e., the efflux from the doublet) 
“pushes” the Chaser away from the obstacle. The resulting potential due to an obstacle 
doublet is 

= - || COS (7i) (6.37) 

IkcoJ 2 
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(a) 3-D Representation of a Doublet Potential 



(b) 2-D Doublet Potential Contours and Gradient 


Figure 6.5. Elementary 2-D Doublet Potential Flow. 


where 7 * is the angle between the Chaser and the obstacle given as 

= tan"^ 

and Ycoi is the relative position of the Chaser with respect to the obstacle. The last 
remaining step is to determine the doublet strength associated with each obstacle. 

Assuming the spacing between multiple objects is sufficiently large such that there is lit¬ 
tle near-field interaction between two or more doublets, each doublet in the workspace 
can be treated independently [40]. Due to the linearity of Laplace’s Equation and sub¬ 
sequent elementary solutions, the strength of obstacle doublet is the superposition of 
the doublet strengths associated with the individual interactions between the obstacle 
doublet and other elementary solutions present in the workspace. That is, the total 
strength of the doublet is the sum of the individual doublets strengths when consid¬ 
ering all the other elements independently and can be written as 

Ki = Kuf + l^tar (6.38) 

where Ki^uf is the doublet strength due to the interaction of the obstacle doublet 
and the uniform flow while Ki^tar is the doublet strength due to the interaction of the 
obstacle doublet and the Target doublet. By applying the von Neumann boundary con¬ 
dition to each obstacle doublet, the strength of the doublet can be found as a function 
of the radius associated with the obstacle keep-out zone. 
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First, consider the interaction between the obstacle doublet and the uniform flow. Simi¬ 
lar to as before, at the boundary of the circular obstacle, the radial velocity is zero and 
Ycd = Ri, yielding an obstacle doublet strength due to the interaction with the uniform 
flow of 

«i,uf = VreiRi COS (ttg) (6.39) 

Similarly, the obstacle doublet strength due to the interaction with the Target doublet is 




^tar 


cos {jt) 

(l|ro.i|l2-^*)'cos (7i) 


(6.40) 


where Vo^t is the relative position of the obstacle with respect to the Target spacecraft. 
Substituting Equation (6.39) and Equation (6.40) into Equation (6.38), the resultant 
obstacle doublet strength is given as 

Ki = Wrefi?- COS (ttg) Ktar 77 - 2 (6-41) 

(||rit||2-i?i)^cos (7i) 


The resulting potential field given the elementary solutions, docking cone corridor con¬ 
straint, and any obstacles is 

K K 

0tot = ^'ref Ike/II 2 COS (a^) + COS (7*) - ^ ^ COS (7^) (6.42) 

The negative gradient of the total potential field in polar coordinates, V0 (r, 9) = 
is given as follows: 

d(J) ^tsr ^ 

-Vr = = -Vref COS (a„) + ^ — COS ( 74 ) - -—^72 COS (7*) (6.43a) 

dr II ret II 2 r llr^oJk 

-ve = --^ = ^'ref sin (a^) + — sin ( 7 *) -k (t*) (6.43b) 

In order to convert the negative velocity potential from a Polar coordinate system in 
Equation (6.43) to a Cartesian coordinate system, the following transformation is uti¬ 
lized: 


No 


- = Riag)V4>l, (r, 9) + R( 7 t)(r, 0) + R(7. - ^t)(r, 9) (6.44) 


2=1 
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where the rotation R(-) is defined as 


cos (■) — sin(-) 
sin(-) cos(-) 


given the angle the Chaser makes to each elementary solution. 


(6.45) 


6.4.3 CL-RRT*-HPF for Spacecraft Guidance and Control 

The proposed sampling-based guidance method, which is given the prefix CL for closed- 
loop, is predicated on the RRT* algorithm (described in Section 6.4.1) and includes the 
capability to bias the configuration space using HPFs as the underlying Navigation 
function. The subtle nuances to this baseline algorithm is the biasing of the random 
samples from the configuration space with a HPF (in a manner similar to [67]) in addi¬ 
tion to utilizing this algorithm in feedback form — thus achieving a “closed-loop”. 

Feedback action is achieved in a manner similar to the IDVD [52], [72] and MPC [52], 
[61] methods by implementing the control to transfer the system from the current (root) 
node to the next node in the resulting path and resolving the path-planning problem. 
Instead of completely rebuilding the tree every guidance sampling period, the proposed 
closed-loop RRT* (CL-RRT*) guidance methods reuse the previous solution, if one was 
found. In a similar fashion to the Dynamic RRT introduced by Fergurson et al. [171], 
the proposed CL-RRT* checks the previous nodes and trajectories for a collision with 
obstacle space Xobs and removes them from the previous solution. The verified nodes 
are then used to initialize a new tree that is rooted at the current position of the Chaser. 

For spacecraft proximity maneuvering, the “default" choice of a cost functional is to 
minimize propellant usage - implying a minimum-propellant optimal control problem 
(described in Section 2.4.1 and 6.2). As such, a minimum AV cost functional was 
chosen as the metric by which to assign the parent node to the randomly sampled 
configuration. That is, the connection between an existing node, g*, and the candidate 
node, gnew which minimizes the Li—norm of the change in velocity required to go from 
Xg to gnew snd from gnew to the goal configuration qj while satisfying the constraints 

Il'T^reff’new.i V*|ll “L H'T^ref (f’/,new fi,new)||]^ (6.46) 

where r-new.i is the unit vector parallel to the line segment connecting g* to gnew and r/,new 
is the unit vector parallel to the line segment connecting gnew to g/. 

In this proposed guidance strategy, a two-fold approach for the avoidance of a dynamic 
obstacle is utilized. Without fitting the motion of obstacle some predefined models, 
straight-line motion of the obstacle is assumed since the underlying motion is rectilin¬ 
ear. However, since this assumption may not necessarily hold, it is desired to adapt a 
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strategy which does not depend upon the velocity of the obstacle. The implication of 
not using this information is to propagate the position of the obstacle is a sufficiently 
fast guidance sampling rate. 

The premise of this two-point strategy is adapted from [38] where the radius of the 
obstacle is artificially increased. In its original implementation, the virtual radius of the 
obstacle, /?', was chosen to be 


R', = Ri + VaAt (6.47) 

where K is the maximum approaching speed and At is the total time delay associated 
with the system. Since RRT-based methods were originally developed for kinematically- 
controlled and driftless vehicles\ it is necessary to adapt it for a dynamically-controlled 
vehicle. 

For a dynamically-controlled vehicle, such as a spacecraft, the smallest control input in 
terms of impulse is governed by the minimum (acceleration) impulse bit of the thruster 

AV = Clthr^^min (6.48) 

where athr is the acceleration due to the control input (i.e., thruster). Rearranging Equa¬ 
tion (6.48) for the minimum time, At^m, and substituting into Equation (6.47) yields the 
virtual radius as a function of the minimum impulse bit of the Chaser: 

R' = R, + ^ (6.49) 

®thr 

Furthermore, the maximum approach velocity, K, is chosen to be 

Ki' = |vll|v,| (6.50) 

where is the velocity of the Chaser and Vj is the instantaneous velocity of the 
obstacle. It is important to note, defining the approach velocity in the traditional sense 
- as the relative velocity between the Chaser and the obstacle - can lead to tempo¬ 
rary conditions where the relative velocity between the Chaser and obstacle are zero, 
thereby causing Ri. This can result in safety concerns where the Chaser moves 
too close to the obstacle causing a possible violation of the obstacle keep-out zone. 
Defining the approach velocity to be Equation (6.50) implies the approach velocity will 
only be zero when either the Chaser or obstacle is stationary. 

^ Recall a kinematically-controlled vehicle Is one whose control Input Is Its velocity. A driftless ve¬ 
hicle Is a vehicle that "can be stopped Instantaneously by setting the control Input to zero" [42] - that 
Is, “momentum-less systems”. Examples of kinematically-controlled and driftless vehicles Include au¬ 
tonomous vehicles [38], [43], [44] as well as robotic manipulators under kinematic control [45H48]. 
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The second point obstacle avoidance is integrated into the proposed strategy is into 
the generation of the HPF gradient. First, the virtual obstacle radius (defined by Equa¬ 
tion (6.49)) can be used to determine the obstacle doublet strength given by Equa¬ 
tion (6.41). Additionally, the effect of the obstacle’s velocity on the gradient is computed 
as 


No 

V0x,tot = V0a; -|- 

i=l 



(6.51) 


where V0a; is the gradient of the static HPF (in the Cartesian coordinate system). The 
gradient in Equation (6.51) can then be directly applied to bias the randomly sampled 
configuration 

Q'rand, biased = *?rand + ( V0a;^tot) (6.52) 


Lastly, when the algorithm terminates (either by successfully completing the path to the 
goal configuration or exceeding the number of iterations) the “solution node”, denoted 
as qs, is chosen as the node that minimizes the distance to the goal configuration 
qj [44]. The solution path is then determined by traversing the tree from the to the root 
node, go- The resulting commanded control input is the one control input associated 
with transferring the system from go to the next node along the solution path - similar 
to how feedback action is achieved for the IDVD [52], [72] and MPC [52], [61] methods. 
The resulting path is then fed-back and the process repeated. 

The remainder of this sub-section will incorporate (dynamic) obstacle avoidance into 
proposed guidance method. 


6.5 Attitude Control 

While each of the algorithms described in Section 6.3 are capable of controlling the 
attitude of the FSS, they are utilized for translation guidance and control only - with 
the exception of the APF and AAPF methods, which also provide attitude guidance 
and control. This was chosen since the single rotational DOF motion of the Chaser is 
largely unconstrained and not coupled to the translational dynamics of the FSS. Two 
forms of attitude control are implemented: a Minimum Fuel Slew plus PD control and 
a saturation controller [172]. A summary of the attitude control methods is listed in 
Table 6.3. 

6.5.1 Saturation Attitude Controller 

Compared to a wide range of controllers, including linear time (in)variant, nonlinear 
time invariance, discontinuous, and optimal controllers, the saturation controller was 
shown by Rao and Bernstein to have a lower maximum under ‘nominal’ conditions 
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Table 6.3. Summary of attitude control methods. 


APF AAPF IDVD LQ-MPC NMPC RRT* RRT*-HPF 
APF-Based x x 

Min. Fuel + PD x x x 

Saturation Controller x x 


[172]. Note, by ‘nominal’, it is meant the conditions under which it was tuned. For 
off-nominal conditions, such as variations in mass, pole perturbations, measurement 
delays, and unmodeled dynamics, the saturation controller was observed to exhibit 
“Good” performance (on a scale of “Good”, “Fair”, “Poor”, and “Very Poor”) [172]. The 
resulting saturation controller implemented is given as. 


M = — sat 


■^^max 


^bJzz^err “h Satg 


Sa'Jzz^err “h I 1 Jzz^e 


(6.53) 


where Sa, Sb e M+; Jzz e M+ is the moment of inertia of the FSS about its vertical axis; 
6*err G M is the error attitude of the FSS; and g M is the error angular velocity of 
the FSS. The saturation controller was tuned (per Reference [172]) with the aid of the 
TRIDENT-GNC nonlinear numerical simulator (described in Section 3.5.1) [31]. The 
parameters chosen for the saturation controller is listed in Table 6.4. 


Table 6.4. Saturation Controller Parameter Selection. 


Parameter 

Value 

Sa 

0.05 s-2 

Sb 

0.0625 s-3 

Jzz 

0.253 kg-m^ 

"rtmax 

0.0437 N-m 


6.6 Experimental Setup 

To conduct the experimental campaign, two FSSs are used - one is the set as the 
Target while the other is the Chaser. The test setup of the POSEIDYN test bed for Case 
1 is illustrated in Figure 6.6 which used a third, third-generation FSS as the obstacle. 
The remainder of the section lists the various guidance parameters used throughout 
the test campaign. 
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Figure 6.6. Test setup for the real-time guidance and control algorithm 
experimental campaign. 

6.6.1 APF and AAPF Parameters 

Relevant guidance parameters are summarized in Table 6.5. Guidance algorithm pa¬ 
rameters were tuned with the aid of TRIDENT-GNC, nonlinear numerical simulator. 
These parameter were then held constant throughout the test campaign in an effort 
to enable greater comparison between the algorithms. It is worthwhile to comment on 
the selection of some of the parameters. First, the potential field shaping matrices, 
Qa, Qd,and P^, were selected to produce minimal chattering upon encountering an 
obstacle boundary constraint using the APF method. Next, the AAPF reference trans¬ 
lational speed, Uref, was chosen to produce a similar rendezvous time compared to the 
Case 0 APF experimental results. The AAPF reference rotational angular speed, Wref, 
was chosen in order to ensure proper alignment prior to entering the docking cone cor¬ 
ridor. Lastly, the adaptive parameter projection limits, p~,pt, were chosen through the 
use of a nonlinear numerical simulator to reduce chattering about the target configura¬ 
tion [31]. 


6.6.2 IDVD Parameters 

Relevant guidance parameters for the IDVD method is presented in Table 6.6 [52], [72]. 

6.6.3 LQ-MPC and NMPC Parameters 

Relevant guidance parameters for the LQ-MPC and NMPC guidance methods are pre¬ 
sented in Tables 6.7 and 6.8, where P is the solution to the algebraic Riccati equation 
associated with the discrete LQR problem [52], [61], [72]. Since the data presented are 
from two different test campaigns prior to the adoption of the proposed Standard Test 
Framework (Refs. [52], [72] and [61]), different tunings are presented for the LQ-MPC 
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Table 6.5. Summary of APF and AAPF guidance parameters. 


Parameter 

Value 

Number of Runs 

5 

Guidance Sample Time, Tg 

0.1 s 

hi 

0.400 m 

0j\ 

1.121 m 

hi 

0.600 m 

0,2 

1.320 m 

0-3 

1.600 m 

k k 

tvai 

1.0 

Qa 

diag([0.025, 0.025, 0.075]) s’^ 

Qd 

diag([0.0125, 0.0125, 0.075]) s'^ 

Pd 

diag([20, 20, 0]) s'^ 

N, 

IsxS 

(Ti 


A 

0.3 

Vrei 

0.025 m/s 

^ref 

2 7s 

Po 

[0.1581, 0,0,1.1581, 0,0.2739]^ 

Pt, Pi 

±0.75 for / = 1,2,..., 5 
±0.25 for / = 6 


Table 6.6. Summary of IDVD guidance parameters. 


Parameter 

Value 

Number of Runs 

1 

Guidance Sample Time, Tg 

0.2 s 

Translation Polynomial Order, 

5 

Time Polynomial Order, ^ 

2 

Maximum Force, Tmax 

0.15 N 

Maximum Time, tf 

100 s 

Maximum Number of IPOPT Iterations 

15 


method. As previously mentioned, experimental results for the rotating hyperplane 
(rotating hyperplane (RH)) and dual hyperplane (dual hyperplane (DH)) methods will 
be presented. 
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Table 6.7. Summary of LQ-MPC guidance parameters for Case 1. 


Parameter 

Value 

Number of Runs 

1 

Guidance Sample Time, Tg 

5 s 

P 

P 

Q 

diag(10^10^10^10®) 

R 

diag(10^, 10^) 

Maximum Force, Tmax 

0.15 N 

Horizon Length, N 

20 

Max. Number of Iterations 

100 

6.8. Summary of LQ-MPC and NMPC guidance paramel 
2. 

Parameter 

Value 

Number of Runs 

3 

Guidance Sample Time, Tg 

3 s 

P 

P 

Q 

diag(l,l,3 x 103,3 x lO^) 

R 

diag(10^, 10^) 

Maximum Force, T^ax 

0.15 N 

Horizon Length, N 

30 

Max. Number of Iterations 

100 


Table 6.9. Summary of RRT* and RRT*-HPF guidance parameters. 


Parameter 

Value 

Number of Runs 

1 

Guidance Sample Time, Tg 

0.1 s 

No. of Nearest Neighbors, k 

20 

Max. Number of Iterations 

250 

Vre\ 

0.025 m/s 


6.6.4 RRT* and RRT*-HPF Parameters 

Relevant guidance parameters for the RRT* and RRT*-HPF guidance methods is listed 
in Table 6.9. In the proposed closed-loop RRT* guidance method, candidate nodes 
select the best parent node which minimizes a Minimum AV cost function (defined in 
Section 6.4.3). 
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6.7 Experimental Results and Discussion 

6.7.1 Results using the Standard Test Framework 

A summary of the averaged translational control effort and maneuver time is presented 
in Figure 6.7. A detailed summary of the averaged translational (and rotational) control 
effort, maneuver time, and the translational component of the non-dimensional control 
effort, z/i for Test Cases 0-3 is presented in Tables 6.11,6.12, 6.13, and 6.14, respec¬ 
tively. Furthermore, the averaged trajectories, distance to the goal time history, and 
(translational) control effort versus the distance to the goal for Test Cases 0-3 is pre¬ 
sented in Figures 6.9, 6.10, 6.11,6.12, respectively. A comparison of the attitude time 
history for each test case is presented in Figure 6.7. Lastly, all methods evaluated suc¬ 
cessfully avoided the obstacle keep-out zone in addition to meeting the attitude error 
requirement prior to entering the docking cone corridor as illustrated in Figure 6.8. 

For the baseline case. Case 0, the two RRT*-based methods performed in a similar 
manner with the AAPF guidance algorithm. These methods were found to exhibit a 
similar linear increase in the translational control effort, as observed in Figure 6.9b. In 
fact, during any obstacle-free portion, this characteristic was observed to hold across 
all test cases for not only this algorithm but also the IDVD method. In contrast, the 
nominal distance time history for the APF resembles that of an exponential response, 
illustrated in Figure 6.9b as well. This is evident as the Chaser crosses the DCC in 
approximately 70 s and spends the remaining 90 s traversing 0.75 m to dock with the 
Target. It is worthwhile to note, the MPC methods exhibit a similar exponential decrease 
in velocity once they are sufficiently close, but to a lesser extent then the APF method, 
as illustrated in Figure 6.1 Ob and 6.11b. 

The ratios of ui to the optimal z/j" are tabulated in Table 6.10. Additionally, for Case 
1, both the IDVD and closed-loop RRT*-HPF (CL-RRT*-HPF) methods were found 
to produce optimal trajectories (compared to the GPOPS - II solution), as illustrated 
by Figure 6.10a. However, while the LQ-MPC/RH method did not produce an opti¬ 
mal trajectory, it did utilize the smallest amount of control effort amongst the evalu¬ 
ated methods. For Case 2, the RRT*-based methods, NMPC, and LQ-MPC/DH were 
found to produce optimal trajectories. For this case, the LQ-MPC method using the DH 
method developed in [61] achieved the smallest expended control effort. Overall, the 
implemented RT-G&C methods were observed to use approximately 5 x -15x greater 
control effort compared to the optimal (GPOPS - II) solution. 
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Figure 6.7. Comparison of the translational control effort and maneuver 
time. 
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Table 6.10. Ratio of the non-dimensional control effort ui for each evalu¬ 
ated method to the optimal non-dimensional control effort z/^. 



Case 0 

Case 1 

Case 2 

Case 3 

APF 

14.17 

14.72 

21.25 

12.88 

AAPF 

10.78 

9.62 

12.38 

11.18 

IDVD 

— 

10.93 

— 

— 

LQ-MPC/RH 

— 

6.88 

8.20 

— 

LQ-MPC/DH 

— 

— 

4.52 

— 

NMPC 

— 

— 

6.06 

— 

RRT* 

10.20 

9.01 

9.83 

— 

RRT*-HPF 

10.35 

9.30 

10.87 

— 


Table 6.11. Case 0 comparison of guidance algorithm performance met¬ 
rics. 



Avg. Translational 
Control Effort 
(N-s) 

Avg. Rotational 
Control Effort 
(N-m-s) 

Avg. 

Maneuver Time 

(s) 

Avg. Translational 

(xlO-4) 

APF 

7.83 ± 0.31 

0.47 ± 0.026 

163.74 ± 3.01 

4.89 ± 0.26 

AAPF 

5.72 ± 0.27 

0.38 ± 0.015 

157.23 ±1.82 

3.72 ± 0.14 

CL-RRT* 

5.60 

0.39 

162.77 

3.52 

CL-RRT*-HPF 

5.53 

0.38 

158.86 

3.57 


Table 6.12. Case 1 comparison of guidance algorithm performance met¬ 
rics. 



Avg. Translational 
Control Effort 
(N-s) 

Avg. Rotational 
Control Effort 
(N-m-s) 

Avg. 

Maneuver Time 

(s) 

Avg. Translational 

(xlO-^) 

APF 

15.29 ± 0.70 

0.76 ± 0.048 

249.89 ±11.69 

6.27 ± 0.29 

AAPF 

12.15 ± 0.50 

0.81 ± 0.067 

302.96 ± 6.36 

4.10 ± 0.19 

IDVD [52], [72] 

4.35 

0.20 

95.16 

4.67 

LQ-MPC/RH [52], [72] 

4.06 

0.22 

141.56 

2.93 

CL-RRT* 

6.82 

0.41 

182.26 

3.84 

CL-RRT*-HPF 

6.10 

0.39 

157.66 

3.96 
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Table 6.13. Case 2 comparison of guidance algorithm performance met¬ 
rics. 



Avg. Translational 
Control Effort 
(N-s) 

Avg. Rotational 
Control Effort 
(N-m-s) 

Avg. 

Maneuver Time 

(s) 

Avg. Translational 

(xlO-4) 

APF 

13.04 ± 1.17 

0.61 ± 0.011 

182.36 ± 5.49 

7.33 ± 0.83 

AAPF 

10.99 ± 0.47 

0.73 ± 0.056 

178.96 ± 2.49 

4.27 ± 0.23 

LQ-MPC/RH [61] 

6.07 ± 1.03 

0.36 ± 0.013 

220.52 ± 7.60 

2.83 ± 0.55 

LQ-MPC/DH [61] 

3.73 ± 0.10 

0.39 ± 0.044 

244.81 ± 5.62 

1.56 ± 0.033 

NMPC [61] 

4.17 ± 0.63 

0.33 ± 0.038 

207.91 ± 33.76 

2.09 ± 0.49 

CL-RRT* 

4.86 

0.33 

147.13 

3.39 

CL-RRT*-HPF 

5.96 

0.40 

162.78 

3.75 


Table 6.14. Case 3 comparison of guidance algorithm performance met¬ 
rics. 



Avg. Translational 

Avg. Rotational 

Avg. 

Avg. Translational 


Control Effort 

Control Effort 

Maneuver Time 



(N-s) 

(N-m-s) 

(s) 

(xlO-^) 

APF 

9.84 ± 0.43 

0.57 ± 0.058 

204.60 ± 6.09 

4.92 ± 0.23 

AAPF 

7.47 ± 0.39 

0.47 ± 0.06 

178.96 ± 2.49 

4.27 ± 0.27 
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lAttitude Errorl (°) lAttitude EiTorl {°) 



0 0.5 1 1.5 2 2.5 3 3.5 4 


Distance Remaining (m) 

(a) Case 0 



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 



Distance Remaining (m) 


(b) Case 1 



Distance Reamining (m) 


Distance Remaining (m) 


(c) Case 2 


(d) Case 3 


Figure 6.8. Comparison of the attitude error versus distance remaining. 
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(a) Case 0: Attitude Error Time History (a) Case 1: Trajectory 




(b) Case 0: Distance to Goal Time History (b) Case 1; Distance to Goal Time History 




Distance Reamining (m) 


(c) Case 0: Thruster ON-Time Profile 


(c) Case 1: Thruster ON-Time Profile 


Figure 6.9. Case 0 experimental 
results. 


Figure 6.10. Case 1 experimental 
results. 
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(a) Case 2: Trajectory (a) Case 3: Trajectory (with Distance to obstacle inset) 



Time (s) 


(b) Case 2: Distance to Goal Time History 



(b) Case 3: Distance to Goal Time History 



(c) Case 2: Thruster ON-Time Profile 


(c) Case 3: Thruster ON-Time Profile 


Figure 6.11. Case 2 experimental 
results. 


Figure 6.12. Case 3 experimental 
results. 
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6.8 Simulation-Based Case Study Evaluating the Effect 
of Dynamic Obstacles on G&C Methods 

A simulation-based Taguchi Analysis^ was performed in order to assess the effect of 
a dynamic environment on the non-dimensional control effort z /1 for various classes of 
RT-G&C methods. The three algorithms considered — the APF, AAPF, and IDVD guid¬ 
ance methods — span the range of RT-G&C from analytical methods (APF) to hybrid 
(AAPF) and optimization-based (IDVD) methods. To assess the effect of a dynamic en¬ 
vironment, the rotation radius and radius of rotation of the obstacle in Case 3 are varied 
per Table 6.15. A 19 ( 8 ^) Taguchi array^ was used for the two-level, three-factor design 
of experiment (DOE), resulting in nine experiments for each guidance method [138]. 
Each experiment was performed using the TRIDENT-GNC numerical simulator [31] 
(described in Section 3.5.1). The normalized average effect of each factor on the per¬ 
formance criterion z/i is presented in Figure 6.13. As one would expect, the smallest 
value of z/i occurs when both the rotation rate and radius of rotation are small. Note, as 
the these factors go to zero. Case 3 effectively becomes Case 0 with a point-obstacle. 

Table 6.15. Summary of the values for the factors considered in the case 

study. 


Level 

Rotation Rate 
(deg /s) 

Radius of Rotation 
(m) 

1 

2 

0.1 

2 

4 

0.25 

3 

8 

0.5 


As illustrated, each of the three methods (representing each class of RT-G&C algo¬ 
rithms) were observed to exhibit an effect on z/i proportional to the rotation rate of the 
obstacle. Both the analytical and hybrid methods, APF and AAPF respectively, exhib¬ 
ited similar sensitivities while the optimization-based method, IDVD, was observed to 
have a greater sensitivity. Furthermore, unlike the hybrid and optimization-based meth¬ 
ods, the analytical method was found to be insensitive to the radius of rotation of the 
obstacle. For the hybrid method, the radius of rotation was found to exhibit an effect on 
z/i proportional to the radius of rotation. However, for the optimization-based method, 
the effect of the radius of rotation on ui is more complex. This complexity may arise 
from possible interactions between the rotation rate and radius of rotation. Additionally, 
the optimization-based method exhibited a larger range of effect on z/i compared to the 
other analytical and hybrid methods. 

^See apprefappendix:taguchMethod for more information on the Taguchi Analysis. 

^See Table B.1 for the structure of the L9(3^) orthogonal array. 
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Figure 6.13. Average effect of varying the rotation rate and the radius 
of rotation of the obstacle for an analytical (APF), hybrid (AAPF), and an 
optimization-based (IDVD) guidance and control methods. 


It is worthwhile to comment on the effect of the sampling rate of the guidance method 
as well as knowledge of obstacle’s motion. If the motion of the obstacle is known and 
can be predicted (with high confidence), this information can be successfully integrated 
into the path-planning process in addition to a larger sampling period can be used. 
When the motion of obstacle is unknown, or predicted with low confidence, a smaller 
sampling period must be used to compensate for the lack of information. A similar trend 
was found to hold for the MPC-based methods. Flowever. since significant tunning of 
the various MPC parameters (e.g., sampling period, Q, R) was found to be required, 
these methods were not considered for this study. 


6.9 Comparison of Guidance and Control Methods 

In order to fully assess the guidance and control algorithms using the Guidance Com¬ 
parison Metric proposed in Chapter 5, two separate plots are required. The first plot 
compares the non-dimensional control effort ui across each test case, illustrated in Fig¬ 
ure 6.14 and 6.15. The second plot than compares the non-dimensional control effort 
ui to the non-dimensional sampling period 1 ^ 2 , illustrated in Figure 6.16 and 6.17. 

6.9.1 Overall Performance Assessment 

The translational component of the non-dimensional control effort, z/i, for each method 
across each test case is illustrated in Figure 6.14 and 6.15. Note, the optimal control 
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solutions found in Section 6.2 are included for reference purposes. Additionally, the 
numerical values for are provided in Tables 6.11,6.12, 6.13, and 6.14. 
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Figure 6.14. Comparison of the APF and AAPF to other guidance and 
control methods tested using the Standard Test Framework. 
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Figure 6.15. Comparison of guidance and control methods using the Guid¬ 
ance Comparison metric. 
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Of the algorithms evaluated, the APF method was observed to have the largest range 
of values, with its largest value being associated with Case 2. In contrast, the AAPF, 
CL-RRT*, CL-RRT*-HPF, and MPC methods were observed to produce smaller varia¬ 
tions across the test cases for which they were evaluated. Additionally, the ui values 
for the optimal control solutions were found to be relatively consistent across all the 
test cases. However, Cases 1 and 3 are larger by approximately 21.0% and 10.2% 
respectively. 

While the AAPF was observed to have a smaller variation across the four test cases, ui 
was observed to trend in a similar fashion to the APF method. This trend is indicative of 
the sensitivity to the number and geometry of the obstacles for the APF-based methods. 
However, despite this observed characteristic of APF-based methods, the use of an 
adaptive feedback law was observed to mitigate this sensitivity significantly for a small 
penalty in the increased complexity of the AAPF method. Similarly, the ui values for 
both of the RRT-based methods were found to trend in a similar fashion. The increase 
in z/i from the baseline (Case 0) to Case 1 for the CL-RRT* method was observed to be 
approximately 8.7% while the CL-RRT*-HPF method was observed to have a 10.4% 
increase. While this increase was smaller than that of the optimal control solution, 
the value of ui is roughly an order of magnitude larger. Compared to the APF-based 
methods, the value of ui was observed to decrease between Case 1 and Case 2 rather 
than increasing. These observations regarding the trending of the APF-based and 
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Figure 6.16. Comparison of guidance and control methods using the Guid¬ 
ance Comparison metric. 
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Figure 6.17. Comparison of guidance and control methods using the Guid¬ 
ance Comparison metric. 


RRT-based methods support the claim regarding the ability of this metric to characterize 
guidance and control methods. 

Lastly, the proposed Guidance Comparison metric, z = is illustrated in Fig¬ 

ure 6.16 and 6.17. The sensitivity of each method to the number and geometry of 
obstacles can be observed by the tightness of the grouping it forms. That is, the APF 
method was observed to generate a loose grouping (i.e. a larger range) of values for 
z/i - indicating a larger sensitivity to number and location of the obstacles - while the 
AAPF and RRT-based methods were observed to form a tight grouping. Additionally, 
while MPC-based methods were observed to attain better performance (i.e. a smaller 
z/i value) than the AAPF and RRT-based methods, it can be observed that both of 
these methods have the ability to react to a more dynamic environment (due to their 
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lower computational time) while maintaining a comparable level of performance to the 
MPC-based methods. 


6.9.2 Assessment of Path-Planning Capability 

From these test cases, one can draw conclusions regarding the path-planning capabil¬ 
ities of each algorithm. Despite the pathological nature of Case 1, both the APF and 
AAPF methods were able to reach the goal position despite their susceptibility to local 
minima. For navigating multiple obstacles, the AAPF method was observed to produce 
a smoother, more efficient trajectory, compared to the APF method. This was attributed 
to a slower changing gradient in the vicinity of the obstacles, as illustrated by the slowly 
decreasing average distance to the goal time history in Figure 6.11 b. Both the APF and 
AAPF methods were observed to successfully avoid the moving obstacle, as illustrated 
by the inset in Figure 6.12a. Note, these methods encountered the obstacle at differ¬ 
ent times and positions. The average closest approach of the Chaser to the moving 
obstacle using the APF method was 0.424 meters and occurred at a distance of just 
over 2 meters from the goal - corresponding to an average obstacle angular position of 
approximately -64° from the initial angular position. Likewise, the average closest ap¬ 
proach using the AAPF method was 0.446 meters occurred at a distance of just under 
2 meters from the goal; however, due to the speed difference of the Chaser, the aver¬ 
age angular position of the obstacle was approximately -158° from the initial angular 
position. Flowever, despite the position differences of the obstacle when it was encoun¬ 
tered, the AAPF method was able to generate a more efficient trajectory compared to 
the APF method. 



X(m) 


(a) CL-RRT* Initial Tree for Case 1 



X{m) 


(b) CL-RRT*-HPF Initial Tree for Case 1 


Figure 6.18. Comparison of the Initial Trees Generated for Case 1. 
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For the RRT*-based methods, the effect of biasing the randomly sampled configuration 
space can be observed by comparing the initial trees for Case 1 in Figure 6.18. The 
primary effect of biasing results in increased directionality of the tree, reducing the 
exploration tendency associated with the construction of the tree [160]. Compared to 
the observed final trajectories, the CL-RRT*-HPF was able to generate a path similar 
to the I DVD method. 

Lastly, a large difference in the path-planning capability of the MFC methods was ob¬ 
served for Case 2, as illustrated by Figure 6.11a. Compared to the LQ-MPC method 
using a rotating hyperplane (LQ-MPC/RH) convexification method, both the LQ-MPC 
using the dual hyperplane (LQ-MPC/DH) method as well as the NMPC method were 
able to resolve the straight-line path. In fact, the LQ-MPC/DH was observed to use 
11.14% less control effort than the NMPC method. 


6.9.3 Effect of the Modulation Technique on ui and z /2 

As illustrated in Chapter 4, the use of the SAM provided superior signal tracking ca¬ 
pability compared to the PWM. This has important implications on the guidance and 
control method and z/i and z/ 2 . First, the SAM provides more accurate realization of 
the commanded input as a result of its feedback loop. Consequently, a slight higher 
control effort is expended - when compared to the PWM - resulting in a larger value of 
z/i. As a direct result of this more accurate command input realization by the SAM, a 
larger guidance sampling period (and larger z /2 value) can possibly used. This implies 
that more computationally intensive algorithms could be used due to the availability of 
a larger sampling period. 


6.10 Summary 

Several guidance and control methodologies were experimentally evaluated using the 
proposed Standard Test Framework for a rendezvous and docking scenario. The op¬ 
timal solution to maneuvering a spacecraft in a cluttered and constrained environment 
was found and presented for each of the test cases associated with this framework. 
Next, several algorithms representing the state of the art, including Artificial Poten¬ 
tial Functions (APFs), Adaptive Artificial Potential Functions (AAPFs), Linear-Quadratic 
(LQ-MPC) and Nonlinear Model Predictive Control (NMPC), as well as Inverse Dynam¬ 
ics in the Virtual Domain (IDVD), were experimentally evaluated and compared to the 
proposed Closed-Loop RRT*-HPF method. The resulting comparison was facilitated 
through the use of the proposed Guidance Comparison metric. Compared to the cur¬ 
rent state of the art, the proposed RRT-based guidance method was found to produce 
near-optimal trajectories and was found to be comparable, in terms of performance, to 
the current state of the art. Furthermore, several key characteristics were demonstrated 
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for the Guidance Comparison metric, including its ability to characterize guidance and 
control methods in addition to its efficacy in comparing guidance and control methods. 
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CHAPTER 7: 

Development of the Guidance and Control Trade 

Space 


The greatest results are simple and may even look trivial 

—Dennis S. Bernstein 
“A student’s guide to research” [173] 


A journal version of the work presented on the proposed Guidance Comparison Met¬ 
ric is under review for publication to the AIAA Journal of Guidance, Control, and Dy¬ 
namics [75]. The work presented in Section 7.3.1 was originally presented at the 27th 
AAS/AIAA Spaceflight Mechanics Conference in San Antonio, TX (February 2017) [74]. 


7.1 Case Study Overview 

To illustrate the extensibility of the proposed characterization framework, this chapter 
presents a case study illustrating the development of a never-before-seen guidance 
and control trade space to aid in the selection of a G&C method for a notional mission. 
Knowledge of the impact of the guidance and control method on the system can inform 
the design earlier in the design process. To aid in this selection, a new and unique 
perspective on the age-old problem of selecting the minimum sampling rate for real-time 
guidance and control algorithms is proposed. Furthermore, several bounding methods 
are presented in an effort to define the boundaries of the trade-space. 

The parameters for this example case study are given in Table 7.1. Note, the reference 
orbit is that of the ISS [77]; the spacecraft parameters are based on a FSS [31], [100]; 
and a typical driving mission requirement, the AV budget, is arbitrarily chosen to be 
within the capability of the system. 


7.2 Bounding the Control Effort 

A driving requirement for autonomous RPO missions is the AV budget. Specifically, 
as the efficiency of the guidance increases, the useful life of the spacecraft is larger 
compared to a less efficient algorithm. Using the parameters in Table 7.1, an upper 
bound on can be estimated given one additional parameter, the maximum accumu¬ 
lated thruster ON-time, AU- Note, in the case that this parameter is not available, it can 
be estimated given these parameters. 
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Table 7.1. Summary of parameters for the example case study. 


Parameter 

Value 

Reference Orbit 

400 km 

(Wet) Mass, mo 

10 kg 

Maximum Thruster Force, T^ax 

0.20 N 

Specific Impulse, Igp 

20 s 

Propellant Mass, mp 

0.5 kg 

At7required 

2.5 m/s 

Executive Time, Tb 

0.1 s 


The change in velocity associated with the linear momentum imparted by a thruster is 
given as 

T M 

(7.1) 

m 

Given that for a spacecraft, the associated control effort considered in the defining z/i 
is the linear momentum impulse, Equation (7.1) can be rewritten as 


AV^ = goMmVi 


(7.2) 


where go is the gravitational acceleration constant and is the maneuvering time. 
For the purposes of upper-bounding z/i, it can be assumed At^ = A^. This assump¬ 
tion allows z/i to be redefined in terms of the required AV and maximum accumulated 
thruster ON-time: 


1^1,u 


At7required 

^hgo 


(7.3) 


Lastly, if Atb is not known, it can be approximately by manipulating the equation defining 
the specific impulse 


Atb 


^p9o^sp 




(7.4) 


Substituting Equation (7.3) into Equation (7.3) yields the upper bound for i>i in terms 
system parameters: 


^lu = 


TtnaxA Vrequired 


rripglhp 


(7.5) 


7.3 Bounding the Sampling Time 

This section will detail several considerations in an effort to bound the sampling time 
and aid in defining the trade space. Section 7.3.1 details a new method on choosing a 
sampling period for a guidance or control task based on first-principles. Section 7.3.2 
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discusses an approach on specifying a time constant associated with the motion of a 
resident space object. 


7.3.1 Minimum Sampling Rate for Real-Time Guidance and Control 


Metric Formulation 

The work in this sub-section was originally presented at the 27th AAS/AIAA Spaceflight 
Mechanics Conference in San Antonio, TX (February 2017) [73]. Modern controllers 
are designed and analyzed in the continuous-time domain and implemented, however, 
in a discrete manner via a digital computer, typically at a fixed interval, Tc. The use 
of a digital computer fundamentally implies a non-zero delay to computational delay. 
Ate, associated with sensor acquisition and computation of the control input. This 
further implies the control output is held constant over the next computational period. 
Furthermore, actuators typically realize the control input through the conversion of a 
digital commanded input to an analog signal via a digital-to-analog conversion. As a 
result of this quantization process, the output of the actuator is limited to discrete levels 
between zero to it maximum output. Exploiting these inherent features of a physical 
system, a minimum sampling rate metric can be derived from first principles. 

Consider first the generic model for a nonlinear system 

x{t) = f{x{t),u{t),t) (7.6) 


Over short the period over which the actuator output remains constant, the motion of 
a mechanical system is assumed to be sufficiently described as rectilinear motion (i.e., 
double integrator dynamics). As a result, the time-averaged acceleration can be written 
as 


X (t + Ate) — x(t) 


^t+Atc 


Ate 


u{t) dr 


(7.7) 


where u{t) is the control input. Note, the first order difference equation in Equation (7.7) 
approaches the derivative / dt = um the limit as Ate 0. Rewriting Equation (7.7) 
as 

Ax = uAte (7.8) 


implies a minimum acceleration impulse which can be imparted onto the system. 


Before continuing, it is worth noting, in the context of defining real-time for a system, 
the acceleration term can be interpreted as a dynamic bandwidth, Cbw- The dynamic 
bandwidth is defined as 


Cbw 


Minimum Control Resolution 
Appropriate Physical Quantity 


(7.9) 
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Analogous to a system bandwidth, cjbw, the dynamic bandwidth is a system-level metric 
which can be used to provide insight into the responsiveness of the system. However, 
unlike the system bandwidth which is typically considered to be inversely proportional 
to the maximum responsiveness of the system, the dynamic bandwidth captures the 
minimum responsiveness of the system [174], Additionally, while the system bandwidth 
is only formally defined for linear systems, the dynamic bandwidth is defined for both 
linear and nonlinear systems [129], [174], 

Continuing, Equation (7.8) implies a corresponding change in the state of the system 
which can be written in terms of the dynamic bandwidth as 

I 

X {t + Ate) - x{t) = Ax =-(ewAtl (7.10) 

Furthermore, given an acceptable change in state (i.e. error) and dynamic bandwidth for 
a particular application, the maximum sampling time, Trt, which captures the dynamics 
of the system can be estimated as 


_ 1 1 ^ 


(7.11) 


Therefore, the resulting minimum sampling rate, /rt for a control system to be consid¬ 
ered real-time is 


/rt 



(7.12) 


Both the maximum sampling period and minimum sampling rate metrics proposed in 
Equation (7.11) and Equation (7.12) respectively have several implications. First, both 
metrics are bounded by the physical system (i.e. mass, inertial , actuator strength) 
as well as the control system implementation. The proposed metrics provide insights 
into the system which are easily mapped to the physical system. Furthermore, these 
metrics provide the designer with ’knobs’ they can adjust in order to estimate the affects 
of changes on various aspects of the system. 

Validation of this metric is presented in Appendix A in an effort to retain the focus of this 
Chapter. 


Guidance Path Re-Planning Rates 

The applicability of this metric is not only for the low-level controller as demonstrated 
in Appendix A; it can also be applied to the guidance, or path-planning, task as well. 
Since optimal control is predicated on a predictive model, uncertainty in the navigation 
solution, actuator dynamics, and unmodeled disturbances can degrade the quality of 
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the solution over time and require an update to the problem. As a result, this metric can 
be used to estimate the rate at which it is necessary to generate a solution. To do so, 
however, requires the estimated error be bounded by some (error) function, e(t), 

e(t)>|/-/| (7.13) 

where / is the dynamics of the system and / is the estimated dynamics of the system. 
Next, differentiating Equation (7.13) with respect to time yields the error dynamics, e. 
Note, the change in error over the real-time sampling period is 

Ae(t) = e(t)TRT (7.14) 


Substituting Equation (7.14) into Equation (7.11) as the change in position, one can 
solve for the time necessary to replan given the bounded error function e(t) as 


= 

2 Cbw 


(7.15) 


7.3.2 Time Constant of a Resident Space Object 

The challenges associated with specifying a time constant of a RSO, Trso, is that it 
requires knowledge of the RSO and the trajectory relative to the the RSO - which 
may be not known a priori. Specifically, two cases arise for specifying a time constant 
associated with a RSO: force and unforced motion. For the case of unforced motion, 
the relative velocity, is proportional to the mean motion,n, of the RSO [78]. Therefore, 
the time constant for the unforced case is can be taken to be the inverse of the mean 
motion 

Trso,u=- (7.16) 

n 


The time constant associated with the forced motion, Trso,f, can be specified by 


Tj-so.f — 


K, 


rel.max 


(7.17) 


where Vjgi.max is the maximum relative velocity between the Chaser and the RSO, and 
Omax is the (estimated) maximum acceleration of the RSO. However, it is unlikely the 
maximum relative velocity will be known a priori. To provide an estimate, Vjgi.max can 
be chosen to be 1 m/s. This simplification allows the time constant to be inversely 
proportional to the estimated maximum acceleration. If the acceleration is unknown, 
the time constant of the RSO can be bounded by using an acceleration from Table 5.3 
or Figure 5.4. 
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7.4 Analysis and Discussion 

Assuming no a priori information is available regarding the forced motion of a space¬ 
craft, the parameters listed in Table 7.1, and a soft-constraint on the acceptable control 
error of 2 m, the resulting time constants and bounds on z/i and z /2 are tabulated in Ta¬ 
ble 7.2. The bounds on the acceleration were chosen to be the maximum acceleration, 
ttmax, and the acceleration given by abar 


r^bar riavg Oj\a 

where aavg is the average acceleration and ai^ is the la standard deviation listed in 
Table 5.3. The time constants specified in Section 7.3 are normalized by dividing by 
the Executive time (Tf,). The resulting bounds are overlaid on the Guidance Comparison 
plot, illustrated in Figure 7.1 to create a guidance and control trade space. Note, the 
grayed-out areas are the infeasible areas based on the bounds on vi and z/ 2 . 

Table 7.2. Summary of time constants and bounds on vi and z/ 2 . 


Parameter 

Value 

Trt 

3.16s 

T'rsO.u 

909 s 

^bar 

0.0164 m/s^ 

^max 

0.0431 m/s^ 

^RSO.f 

{ 23 . 2 , 61 } s 


5.20 X 10 “^ 

1^2 

{ 31 . 6 , 250 , 9090 , 12340 } 


As illustrated, the upper bound on z/i eliminates the APF method as a possible guid¬ 
ance choice due to its high averaged control effort required. For the RSO time constant 
associated with a^ax, the relatively large acceleration of Hayabusa eliminates the MFC 
methods as a possible guidance choice, as illustrated in Figure 7.1 a. If the acceleration 
of the RSO is taken to be abar, Trso,u becomes larger and the MFC methods become 
a feasible guidance method illustrated in Figure 7.1b. Regardless of the bounds cho¬ 
sen for Trso,u, the RRT-based methods and AAPF methods - and possibly the IDVD 
method - are feasible methods for this particular application. Given these results, it 
is clear the RRT-based methods yield a slightly smaller value for z/i compared to the 
AAPF methods. However, it is worthwhile to note, this methodology does not take into 
account any implementation complexities - which is more or less a qualitative assess¬ 
ment - associated with the -potential methods. 
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Figure 7.1. Annotated guidance and control trade space based on the 
Guidance Comparison metric. 
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CHAPTER 8: 

Summary and Future Work 


When we are engaged in research, especially when we are 
working on our dissertations, we think that research is everything, 
but it is not. There is joy in discovering a new result, but I think 
research (in Engineering especially) is most satisfying when it 
serves some immediate practical purpose as well. 

Research should be a source of joy, of exhilaration, and, in many 
ways, an act of love. If it isn’t, then it may be difficult to endure the 
hardships that research entails. 

—Malcolm D. Shuster 
“Advice to Young Researchers” [175] 


8.1 Summary of Contributions 

The work presented in this dissertation focuses on the advancement of the state of the 
art in the evaluation of Real-Time Guidance & Control (RT-G&C) and Computational 
Guidance & Control (CG&C) methods. This work is fundamental to the emerging field of 
CG&C algorithms as it provides the necessary methodology and tools to evaluate their 
performance in an equitable manner. This was accomplished by: significantly improv¬ 
ing the TRIDENT-GNC (Toolkit for Real-time DevelopmENT of Guidance, Navigation 
& Control) software architecture for the Naval Postgraduate School (NPS) POSEIDYN 
(Proximity Operation of Spacecraft: Experimental hardware-ln-the-loop DYNamic simu¬ 
lator) test bed; performing an extensive characterization of the POSEIDYN test bed; an 
experimental evaluation framework and benchmarking test scenarios were established; 
and a global guidance comparison metric was developed and validated. Additionally, 
the challenge of efficiently modulating the guidance control signal into a sequence of 
discrete, binary signals which control spacecraft thrusters was explored - effectively 
solving the low-level actuator control problem. The fully characterized test bed was 
then utilized to compare a newly proposed RT-G&C algorithm based on the Rapidly- 
Exploring Random Tree Star (RRT*) method to existing algorithms in literature using 
the evaluation framework. Lastly, the global Guidance Comparison Metric was then 
utilized to compare the evaluated algorithms. 
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8.1.1 Extensive Characterization of the NPS POSEIDYN Testbed 

A detailed description of the Naval Postgraduate POSEIDYN test bed has been pro¬ 
vided. The quasi-frictionless and low residual acceleration dynamic behavior of the 
test vehicles operating on the POSEIDYN test bed are limited to two translational and 
one rotational degree-of-freedom motion. Multiple vehicles can be operated simultane¬ 
ously on the POSEIDYN test bed allowing research on coordinated control of spacecraft 
teams. The on-board actuators, composed of eight thrusters and/or a reaction wheel, 
are equivalent to the actuators found in orbital spacecraft and further enhance the dy¬ 
namic equivalence. With an available simulator and extensive software development 
tools, guidance and control algorithms for real-time execution on-board the test vehi¬ 
cles can be quickly and easily developed. From the characterization performed of the 
POSEIDYN test bed, the real-time operating system was shown to have bounded laten¬ 
cies and low-residual acceleration of 19 /ig — a claim no other test bed can make. Used 
as a last stage of on-the-ground validation prior to on-orbit deployment or simply as a 
more realistic development environment for novel guidance and control approaches, 
this state-of-the-art dynamic hardware-in-the-loop test bed will continue to be fruitful for 
the advancement of spacecraft proximity operations research. 

8.1.2 Development of the TRIDENT-GNC Software Architecture 

The TRIDENT-GNC (Toolkit for Real-time DevelopmENT of Guidance, Navigation & 
Control) software architecture is composed of a (numerical) development simulator of 
the floating spacecraft simulator (FSS), a software template, and custom Simulink li¬ 
brary. This library contains common software, such as the navigation and control sub¬ 
systems which is used in both the simulator and the FSS autogenerated onboard soft¬ 
ware. This commonality between the simulator and FSS onboard software allows for 
rapid development of multi-rate Guidance, Navigation, and Control (GNC) algorithms in 
a simulation environment and software generation for use onboard the FSS for testing. 

8.1.3 Development of a Standardized Test Framework to Evaluate 
Autonomous G&C Algorithms 

The proposed Standard Test Framework defines relevant scenarios that test various ca¬ 
pabilities of candidate algorithms in a relevant and dynamic environment. This frame¬ 
work is composed of four test cases. The first case is a straight-in approach with 
only the docking cone corridor (DCC) and aids in establishing the baseline. The next 
two cases test the path planning capabilities of the candidate algorithm with a static 
obstacle. The last two cases test the ability of the algorithm to adapt to a dynamic 
environment with one and three obstacles rotating about some predefined point in the 
workspace. The Standard Test Framework is used in Chapter 6 to evaluate and com¬ 
pare several guidance and control algorithms. 
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8.1.4 Development and Validation of a Global Guidance Compari¬ 
son Metric 

The proposed global Guidance Comparison Metric provides a unique capability to com¬ 
pare different algorithms across different architectures. This proposed metric, never 
seen before in the literature, allows for an unbiased comparison of algorithms across 
multiple system architectures through the use of a vector-valued metric, z = [z/i,z/ 2 ]^. 
Due to its generic formulation, the Guidance Comparison Metric can be extended to vir¬ 
tually any mechanical system. A case study using the Guidance Comparison Metric to 
develop a never-before-seen guidance and control trade space illustrated in Chapter 7. 

8.1.5 Experimental Evaluation of the Sigma-Delta Modulator as a 
Thruster Modulation Technique 

Motivated by the inherent properties of a Sigma-Delta Modulator (SAM), an experimen¬ 
tal campaign was conducted to examine its viability as a spacecraft thruster modulation 
technique. The experimental campaign consisted of a rest-to-rest, 90 deg slew ma¬ 
neuver using the reaction control jets onboard a floating spacecraft simulator. In this 
campaign, the SAM was compared against a traditional spacecraft thruster modula¬ 
tion technique, the Pulse-Width Modulator (PWM). The SAM was found to achieve 
a smaller steady-state error in a shorter amount of time, while using less propellant. 
However, due to the short period of the limit cycle induced by the feedback in the SAM, 
the propellant consumption during steady was observed to linearly increase at a higher 
rate than the PWM. Given a set of mission requirements, the SAM can be coupled 
with a deadzone in order to increase the period of the limit cycles, thereby reducing 
the fuel consumption. As an alternative to a deadband, one could alternatively uti¬ 
lize the SAM throughout the transient of a response and switch to a lower-bandwidth 
modulator, such as the PWM, when in the vicinity of the reference state. Additionally, 
candidate closed-loop transfer functions describing the behavior of spacecraft attitude 
control system with a modulator-in-the-loop was developed for both modulation tech¬ 
niques. The candidate SAM model was found to exhibit similar characteristics to both 
the experimental and a simulated nonlinear response. Resultantly, this model, which is 
considered to be validated, can be used in conjunction with linear analysis techniques 
to estimate the response of system. 

8.1.6 Application of a RRT-Based Algorithm as a Feedback Con¬ 
troller for the G&C of a Spacecraft in a Dynamic Environment 

A new guidance and control strategy leveraging the path-generating capabilities of 
Rapid-exploring Random Trees (RRT) as well as the obstacle avoidance capabilities 
of Harmonic Potential Functions (HPF) has been developed and experimentally tested 
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in a rendezvous and docking scenario. To accomplish this, the potential field was de¬ 
veloped using elementary solutions to the Laplace Equation which have strong ties to 
incompressible and irrotational flow. The strengths of the various elements composing 
the potential field were then autonomously assigned given the capabilities of the space¬ 
craft. Next, the gradient of the resulting fields was used to bias a randomly sampled 
state in the workspace towards the goal while remaining in the obstacle-free portion of 
the workspace. Lastly, the path-generating capabilities were utilized to develop a path 
which minimizes the fuel consumption of the spacecraft while reducing computational 
complexity compared to traditional methods to solve optimal control problems. The 
Naval Postgraduate School (NPS) POSEIDYN was then used to validate the real-time 
capabilities of the proposed guidance and control methodology in a computationally 
constrained environment and assess its robustness to realistic noise, delays and un¬ 
certainties. While the two methods provided similar results the straight-path solution, 
the advantage of biasing the random samples from the configuration space is illustrated 
in the single static obstacle case where the chaser must circumnavigate an obstacle. 

8.1.7 Development and Validation of a Metric to Estimate the Min¬ 
imum Sampling Rate 

When implementing a controller for a physical system, the rate at which the system is 
sampled and controlled must be determined. Often times, this rate is chosen through 
various “Rules of Thumbs” which are rooted in the frequency domain. However, a new 
approach to determine the minimum sampling rate rooted in first principles for the real¬ 
time control of a spacecraft was presented and validated via numerical case studies. 
This approach exploits the discrete nature associated with the digital implementation 
of a modern control system along with a unique view on the underlying dynamics to 
derive an estimate for the minimum sampling rate necessary for a control system to 
be considered real-time. This metric, unlike other methodologies, is valid for both lin¬ 
ear and nonlinear systems. Additionally, the resulting formulation easily maps to the 
physical system. Several case studies were performed using relevant linear and non¬ 
linear systems including a spacecraft undergoing a 3-axis reorientation maneuver us¬ 
ing an Eigenaxis-Quaternion feedback controller to a double integrator system under 
a minimum-time control policy. Each system considered utilized a realistic actuator 
with discrete output levels. From these case studies, the proposed metric was found 
to sufficiently estimate the average steady-state output of the system. Additionally, by 
exercising the flexibility of the definition dynamic bandwidth, an upper bound on the 
minimum sampling rate can be found by estimating transient tracking capabilities of the 
system. Furthermore, an extension of this metric was presented to estimate the rate at 
which the guidance path re-planning task must be performed. This metric, while does 
not exactly predict the necessary sampling rate, provides an estimate to well within an 
order of magnitude with few assumptions under ideal conditions. 
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8.1.8 Development and Formulation of a Guidance and Control 
Trade Space for Spacecraft Design and GNC Selection 

Using the experimental evaluation framework and the Guidance Comparison Metric, 
z = a never-before-seen guidance and control trade space was formulated. 

Knowledge of the impact of the guidance and control method on the system can inform 
the design earlier in the design process. This trade space was developed using param¬ 
eters typically available early in the design process. Using these parameters, bounds 
were placed on the non-dimensional control effort., z/i, and non-dimensional sampling 
period, z/ 2 , to narrow down the trade space to select a candidate algorithm. 


8.2 Contribution Traceability 

The supporting evidence for each contribution described in Section 8.1 is cross-referenced 
with the corresponding chapter in Table 8.1. 


Table 8.1. Mapping of the contribution to the supporting chapters. 


Contribution 


Dissertation Chapter 



3 

4 

5 

6 

7 

Extensive characterization 
of the NPS POSEIDYN testbed 

/ 

-- 

-- 

-- 

-- 

Development and implementation 

of the TRIDENT-GNC software architecture 

/ 

-- 

-- 

-- 

-- 

Development of a standardized test framework 
to evaluate autonomous G&C algorithms 

-- 

-- 

/ 

-- 

-- 

Development and validation of a 
global Guidance Comparison metric 

-- 

-- 

/ 

/ 

/ 

Experimental evaluation of the Sigma-Delta 
modulator as a thruster modulation technique 

-- 

/ 

-- 

-- 

-- 

Application of a RRT-based algorithm as a 
feedback controller for the G&C of a 
spacecraft in a dynamic environment 

-- 

-- 

-- 

/ 

/ 

Development and validation of a metric to 
estimate the minimum sampling rate 

-- 

-- 

-- 

-- 

/ 

Development and formulation of a G&C trade 
space for spacecraft design and GNC selection 

-- 

-- 

/ 

/ 

/ 


8.3 Future Work 

In this section, recommendations for future work is suggested. 
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8.3.1 Simulation-Based Standard Test Framework 

The Standard Test Framework proposed in this dissertation provides a framework for 
the experimental evaluation of RT-G&C algorithms. However, test facilities such as 
the POSEIDYN testbed involve a large amount of resources in terms time and capital. 
To make the evaluation methodology described in this dissertation more accessible, a 
simulation-based framework needs to be developed with relevant scenarios and con¬ 
straints. The resulting framework and evaluation methodology can then be adapted for 
more widespread use. 


8.3.2 Spacecraft Control 

While the Sigma-Delta Modulator (SAM) was shown to achieve a smaller steady-state 
error while using less propellant when compared to the Pulse-Width Modulator (PWM), 
the propellant consumption after achieving steady-state was observed to increased at a 
higher rate due to the shorter period of the limit cycle. Future work can explore various 
techniques to reduce the propellant consumption during steady-state while maintaining 
smaller tracking errors (compared to the PWM). Potential avenues of research include 
the incorporation of a deadzone to the SAM, switching to a lower-bandwidth modulator, 
and even exploring more advanced modulators, such as an adaptive SAM [176]. 

8.3.3 Evaluation of G&C Methods 

In addition to evaluating different G&C methods to create a larger repository of data 
points, the G&C methods should be run at various sampling periods. While the RRT*- 
based was shown to produce a near-optimal path, the normalized control effort, ui, 
was observed to achieve a similar value as the AAPF. Therefore, it would be bene¬ 
ficial to determine if there exists any dependency on the sampling period and z/i. If 
a dependency is found to exist, then an optimal value for the sampling period can be 
found which will further reduce the fuel consumption while maintaining a desired level 
of responsiveness. 

8.3.4 Increasing the Computational Performance of G&C Methods 

The ideal optimal G&C method solves the constrained optimal control problem (OCP) 
sufficiently fast. “Sufficiently fast” can be achieved through multiple means including 
computational architectures, speed of the compute element(s), efficient solvers, or bet¬ 
ter predictive models. The first method which should be explored is the use of cus¬ 
tom programming solvers to create more-efficient solvers, such as the those presented 
in [177]-[179]. The use of a custom solver can result in speedup lOx-IOOOx by taking 
advantage of the specific structure of the problem [177]. 
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The second avenue of research to increase the computational performance of G&C 
methods is the use of new compute architectures. Current compute architectures typi¬ 
cally consist of a single serial processor that is shared with other processes. One possi¬ 
ble compute architecture (loosely based on the NVIDIA Tegra series of processors and 
embedded systems, such as the NVIDIA Jetson TK1 [105]), illustrated in Figure 8.1, 
includes the use of a multi-core CPU, co-processor (or secondary CPU), and a graphic 
processing unit (GPU). The multi-core CPU could for a multi-threaded flight software 
design where the GNC algorithms and other housekeeping software run on individual 
threads. 

Lastly, the third avenue of research is the use of better predictive models that can be 
incorporated into the guidance method. This area of research is considered to be an 
extension of the Navigation subset of algorithms of the GNC architecture. It is be¬ 
lieved this can be achieved by using an Interacting Multiple Model filter approach [180]. 
The resulting model probabilities can then be used to inform the guidance algorithm 
of the appropriate model to use. However, the challenge arises of performing this mo¬ 
tion association in an efficient manner as the Interacting Multiple Model approach is 
a resource-intensive task (in terms of computational processing and memory require¬ 
ments) as there n number of filters being run updated at each measurement update. 
The resulting model probabilities can then be used to inform the guidance algorithm of 
the appropriate model to use. 

8.3.5 Extension of the CL-RRT* to RPO 

From the work performed in this dissertation, the advantages of sampling-based meth¬ 
ods are clearly illustrated. Given the higher-dimensionality and nonlinear dynamics 
associated RPO, it would be advantageous to formally extend this method to this do¬ 
main. While this may seem like a daunting task, only two functions in the algorithm 
need to be updated: the steering function which connects candidate nodes to existing 
nodes in the tree and the collision detection routine. 
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Figure 8.1. Notional computational architecture consisting of a CPU, Co¬ 
processor, and GPU. 
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APPENDIX A: 

Validation of the Minimum Sampling Rate Metric 


The work presented in this appendix was originally presented at the 27th AAS/AIAA 
Spaceflight Mechanics Conference in San Antonio, TX (February 2017) [74], 


A.1 Numerical Case Study Overview 

To aid in the validation of this metric, several numerical case studies are performed us¬ 
ing representative mechanical systems, controllers, and actuators with discrete output 
levels. The systems considered include a three-axis reorientation maneuver of a space¬ 
craft using reaction wheels and an Eigenaxis-Quaternion feedback controller [82]; a 
simple harmonic oscillator with a PD controller [129], [174]; and a double integrator (Dl) 
system with reaction control jets using the minimum-time optimal control policy [76]. 

Unless otherwise specified, the mass of the system was chosen to be unity for sim¬ 
plicity. Each actuator is assumed to have a maximum control output of of unity and a 
resolution of 0.001, except for the reaction jets which is a purely discrete actuator with 
two states: ON or OFF. Additionally, each system was tuned to settle within 60 sec¬ 
onds and achieve stability with a delay of 1 second. Furthermore, each system utilizes 
continuous dynamics propagated using a variable step solver and a controller which is 
sampled with a fixed period, T^. The controller sampling period was varied between 
0.001 seconds (1 kHz) to 1 second (1 Hz). Each scenario was run for 120 seconds and 
the average steady-state error in x was recorded over the last 60 seconds. The results 
of each case study were then compared with the predicted minimum change in state 
given by Equation (7.10). To further aid in comparison, the dynamic bandwidth in Equa¬ 
tion (7.10) is replaced with the second-order equivalent system bandwidth, Wbw = ^/m 
and the result is presented. 


A.1.1 Overall System Responses 

The combined responses of each system are illustrated together in Figure A.1. Each 
system was found to exhibit a log-linear response in the average minimum change in 
the state x compared to the sampling rate. Additionally, in order to reduce the numerical 
noise and increase clarity, each response was curve-fit to a log-linear response. This 
response, and its estimated la error are illustrated in subsequent figures to further 
highlight the trends of various responses. 
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Figure A.1. Average steady-state error versus controller sample rate. 

A.1.2 Eigenaxis-Quaternion Feedback Controller 

The first case considered is a three-axis reorientation maneuver of spacecraft using 
an Eigenaxis-Quaternion feedback controller. Additionally, to produce the torque, a 
set of 3 orthogonal reaction wheels were modeled and placed in the loop This system 
is relevant as it is nonlinear in both the system dynamics and control. The system 
dynamics are governed by Euler’s Equation and given as [82], 

T = Jci; + Jlu (A.1) 

where r is the control torque, J is the inertia matrix about the center of mass in the 
spacecraft body Cartesian coordinate system, uj is the angular velocity of the space¬ 
craft body frame with respect to the inertial frame, and is the metrical representa¬ 
tion of the cross product operator [82]. Additionally, the Eigenaxis-Quaternion feedback 
controller is given as 


u(t) = -2kr]e„Jee„ - cJcUgrr + OJgrrJ^err (A.2) 

where Sgrror e is a column vector containing the vector components of the error 
quaternion , r]e„or G M is the scalar component of the error quaternion, and the scalar 
gains c,k e M+. The stability of this controller was proven via Lyanpunov’s Direct 
method [82], [90]. Lastly, the error angle was defined as [82]: 

terror ~ 2 COS ("terror) (A.3) 
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In order to examine the effect of various physical properties, the inertia matrix was 
chosen to be J = diag(l, 0.1,10). Additionally, to meet the 60-second settling time, 
the scalar gains were chosen \o be k = 0.3 and c = 1.0. The results from the case 
study are presented in Figure A.2. When compared to the response of the system, the 
dynamic bandwidth closely matches the first-order fitted response of the numerical ex¬ 
periment. Additionally, the dynamic bandwidth was observed to maintain a similar level 
of error as the inertia was varied between the three axes. Lastly, the predicted minimum 
sampling rate associated with using the system bandwidth and dynamic bandwidth both 
over predict the minimum sampling rate in order to achieve a specific error level. The 
prediction associated with the dynamic bandwidth estimate is within an order of magni¬ 
tude as the numerical response compared to using the system bandwidth metric. For 
this particular case, the predicted error for a given sample time agrees very well with 
the numerical experiment. 

A.1.3 Simple Harmonic Oscillator with PD Control 

The second case considered is a simple harmonic oscillator (SHO) with a PD controller. 
The dynamics of a SHO are given as. 


X = - X (A.4) 

m 

where /c G M+ and is typically referred to as the spring constant for a mechanical 
system [129], [174]. This is a representative second-order system with a common 
controller. As illustrated in Figure A.3 the dynamic bandwidth sufficiently estimates 
the sample time required to achieve a desired level of system performance. Similar 
to the first case, both metrics over-predict the sampling rate necessary to achieve a 
given level of steady-state error. Likewise, the prediction associated with the dynamic 
bandwidth is smaller than that associated with the system bandwidth and is within an 
order of magnitude of the average numerical response. 

A.1.4 Minimum-Time Optimal Control for a Double Integrator Sys¬ 
tem 

The third case study considered is a double integrator (Dl) system under a minimum¬ 
time optimal control policy using discrete actuators, such as reaction jets. There are 
two primary methods of implementing this particular control policy. The first method is 
a clock-based method which assumes knowledge of the switching time. The second 
method (the method implemented for this case study) involves the use of a full-state 
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Figure A.2. Average steady-state error versus controller sample rate for a 
3-axis spacecraft reorientation maneuver. 


switching curve to determine the switching point given as [76] 

Tfl 

S {x{t),x{t)) = x{t) + 77rTl^(^)l^(^) (A-5) 

l\u\ 

where m e M+ is the mass and m e M is the control input. As illustrated by Figure A.4, 
the dynamic bandwidth sufficiently estimates the predicted steady-state error given a 
sampling time. Note, since the system bandwidth and dynamic bandwidth are the 
same numerically, both methods produce the same prediction. Lastly, compared to the 
SHO, the Dl system with discrete actuators produces substantially more variation in the 
steady-state error. 
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Figure A.3. Average steady-state error versus controller sample rate for a 
simple harmonic oscillator. 



Figure A.4. Average steady-state error versus controller sample rate for a 
double-integrator system under the minimum-time optimal control policy. 

A.2 Minimum Sampling Rate for Tracking Transient Re¬ 
sponses 

As demonstrated, this metric provides an order of magnitude estimate for the sampling 
rate required to track a control signal to within a specified level of error. The dynamic 
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bandwidth as given in Equation (7.9) provides the lower-bound on the steady-state 
tracking capability of the system. Redefining the dynamic bandwidth to be equivalent 
to the maximum acceleration of the system 

Maximum Control Effort 
Appropriate Physical Quantity 


one can upper bound the transient tracking capability of the system. Rewriting Equa¬ 
tion (7.8) and substituting in Equation (A.6), the maximum sampling time Trt and mini¬ 
mum sampling rate /rt can be upper-bounded respectively as 


Trt — 


/rt 


Cbw 


(A.7) 


For example, consider the maximum change in error, of the previously ex- 

mained SHO system during its transient period consisting of the first 60 seconds. As 
illustrated by Figure A.5, the dynamic bandwidth provides an estimated sampling rate 
to within 36% of the response of the system, while the system bandwidth significantly 
underestimates the required sampling rate to achieve a desired level of tracking error. 



Figure A.5. Maximum change in error during the transient phase for a 
simple harmonic oscillator. 
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APPENDIX B: 
Taguchi Analysis Method 


B.1 Overview 

The Taguchi analysis method (also referred to as the Taguchi parameter design ap¬ 
proach by some authors [140]) was initially designed for use in process control product 
quality improvements [138], [139]. According to Roy, Taguchi’s most significant contri¬ 
bution was the underlying framework for designing experiments and subsequent anal¬ 
ysis of their results [138]. Due to this wide-reaching contribution, this method has been 
adapted to other applications, such as robust design optimization and analysis [137], 
[140], [181] and experiments involving agriculture and pharmaceuticals [138]. 

The Taguchi method, as it applies to analysis, will be described in this appendix. 


B.2 Methodology 

Roy defines a design of experiment (DOE) as the "technique of defining and investi¬ 
gating all possible conditionals in an experiment involving multiple factors" [138]. It is 
important to note, an input to the DOE is referred to as a factor while the value it takes 
on is referred to as its “level” [137], [138] Conventional approaches have utilized either a 
full-factorial or fractional-factorial unstructured DOE where the number of experiments 
(consisting of a combination of factors at various levels) grows exponentially and the 
contributions of each factor cannot be identified [137], [138], [140]. While the Taguchi 
method does not completely solve the dimensionality issue associated with a DOE with 
many factors, it standardizes both the implementation of the DOE and computation and 
interpretation of the results [138]. 

Compared to the brute-force approach of “think, try this, think, try that”, the Taguchi 
method is applied in four distinct steps: brainstorming, experimentation, analysis, and 
confirmation [138]. To apply the Taguchi method, the experiment must be developed 
using a standard set of orthogonal arrays, referred to as Taguchi orthogonal arrays, that 
specify the minimum number of experiments given the number of factors and the levels 
of each factor [137], [138]. If determining interactions between factors is important, 
the standard Taguchi arrays must be modified. Since interactions are not considered 
in the applications presented in this dissertation, the reader is referred to [137]-[139] 
for further information on interactions and how to modify the standard methodology. 
Since the Taguchi method was formed around measuring deviations from a nominal 
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value [138], [139], the desired output from each experiment is a measurement. This 
measurement, denoted here as the performance criterion (Roy refers to it as an “Overall 
Evaluation Criteria” [138]) and is chosen as the variable which is affected by the level 
of the factors [138]. Lastly, the appropriate level of analysis must be determined. 

As the original application of the Taguchi method was for quality engineering, the anal¬ 
ysis associated with the standard Taguchi method allows the user to identify the impact 
of each factor on the performance criterion, but also identify any interactions (if the 
experiment was chosen to include interactions) as well as pick the set of levels for 
each factor which optimize a certain performance criterion [137], [138]. However, for 
this application, it is only necessary to determine the average effect of each factor on 
the output. As such, the reader is encouraged to refer to [137]-[139] for further anal¬ 
ysis using the Taguchi method. Once the factors, levels, the performance criterion 
are selected, the brainstorming phase is completed [138] and the experiment and data 
collection phase can begin. 


B.3 Taguchi Orthogonal Arrays 


The experimentation phase of the Taguchi method is predicated upon the use of a 
Taguchi orthogonal array. Each Taguchi array is denoted using the following nomen¬ 
clature: 


L 


Afe: 


(</” 


) 


where N^xper is the total number of experiments and Ni^i is the number of levels for 
each of the number of factors, Nfac- In this dissertation, a two-level, three-factor, nine- 
experiment array was selected for the studies performed in Section 5.5.1 and 6.8. The 
19 ( 3 ^) standard orthogonal array used throughout this dissertation is given in Table B.1. 
In this array, the first column corresponds to each experiment, or run. The columns in 
the middle array represent the level for each factor, where each factor is denoted using 
a capitalized letter, beginning with “A” for the first factor. The last column is composed of 
the result, or performance criterion, from nine individual experiments, denoted as Xj, for 
i = 1,2,..., 9. When conducting non-virtual (i.e., physical) experiments, Roy suggests 
running the experiments in a random order to remove any effects due to the experiment 
setup [138]. Once the experiments have been performed and the data recorded, the 
average effects of each factor on the performance criterion can be assessed. 
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Table B.1. 19 ( 8 ^) Taguchi orthogonal array. 



Factor A 

Factor B 

X 

1 

1 

1 

Xi 

2 

1 

2 

X 2 

3 

1 

3 

X 3 

4 

2 

1 

X 4 

5 

2 

2 

X 5 

6 

2 

3 

Xe 

7 

3 

1 

X 7 

8 

3 

2 

Xs 

9 

3 

3 

X 9 


B.4 Computing the Average Effects 

In order to assess the effect, or contribution, of a various factor on the performance 
criterion, the average effect can be simply computed. The average effect of each factor 
at each level can be computed by summing all the elements Xi that contain that factor 
at the level [137], [138]. For example, using Table B.1, the average effect of the factor 
A at level 1 can be computed as follows: 


X1 + X2 + X3 

/ii — - 

3 

Likewise, the computation of the average effect of the factor B at level is 


(B.1) 


X1+X4 + X7 

-Dl = - 

^ o 


(B.2) 


Lastly, each of the average effects may be normalized by the grand average X^vg [137], 


^ ^exper 

Xayg = — ^ Xi 

■‘■^exper 


(B.3) 


Once these values are computed, they can be plotted for visual inspection and inter¬ 
pretation using a quality characteristic (“nominal is best”, “bigger is better”, or “smaller 
is better”) [137], [138]. 


The utility of the Taguchi method arises from the fact that vast insight into the underlying 
sensitivities of the system being analyzed can be made despite the small number of 
experiments performed [137]. To illustrate this, an average effect plot for a notional 
DOE performed using a L 9 ( 3 ^) standard orthogonal array is presented in Figure B.1. 
As illustrated, as factor B is varied from Bi to B^, the performance criterion was found to 
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decrease while the opposite holds for factor A as it is varied from Ai to A 3 . In contrast, 
no effect on the performance criterion is observed when factor A is varied between A 2 
and ^3 as the line connecting the two points is horizontal. Furthermore, a qualitative 
assessment on the presence of any interactions can be determined from this plot simply 
by examining the relatively slopes of the lines for each factor. If they placed on top of 
each other, such as in Figure B.2, one can assess that an interaction between the two 
factors does exist since they are not (nearly) parallel [137], [138]. Further analysis can 
be performed to determine the extent of any interaction between the two factors and 
the average effect on the performance criterion [137], [138]. 



Figure B.1. Notional average effect plot for a 19 ( 8 ^) DOE. 
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Figure B.2. Factor A and B average effects overlaid for the notional 19 ( 8 ^) 
DOE. 
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APPENDIX C: 
Publication History 


C.1 Patents 

J. Virgili-Llop, R. Zappulla II, and M. Romano, “Dynamically tilting flat table to impart 
a time-varying gravity-induced acceleration on a floating spacecraft simulator,” U.S. 
Patent Pending 62 455 775, Feb. 07, 2017 


C.2 Supplemental Media 

Demonstrations of select maneuvers are available on the NPS Spacecraft Robotics 
Laboratory YouTube Channel: 

https://www.youtube.eom/channel/UCyPOVZ2Xrcz_onUMs578aCw. 

A MathWorks user story about the POSEIDYN test bed and the TRIDENT-GNC soft¬ 
ware architecture was published on May 19, 2017 and can be found by clicking here. 


C.3 Journal Publications 

C.3.1 Published 

R. Zappulla II, J. Virgili-Llop, H. Park, C. Zagaris, and M. Romano, “Dynamic air¬ 
bearing hardware-in-the-loop testbed to experimentally evaluate autonomous space¬ 
craft proximity maneuvers,” Journal of Spacecraft and Rockets, 2017, DOI: 10.2514/1 .A33769 
(In Press) 

J. Virgili-Llop, C. Zagaris, H. Park, R. Zappulla II, and M. Romano, “Experimental eval¬ 
uation of model predictive control and inverse dynamics control for spacecraft proximity 
and docking maneuvers,” CEAS Space Journal, 2017, (In Press, 22 May 2017) 


C.3.2 Accepted 

R. Zappulla II, J. Virgili-Llop, and M. Romano, “Spacecraft thruster control via sigma- 
delta modulation,” Journal of Guidance, Control, and Dynamics, 2017, (Accepted 13 
May 2017) 
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C.3.3 Under Review 

R. Zappulla II, H. Park, J. Virgili-Llop, and M. Romano, “Real-time autonomous space¬ 
craft proximity maneuvers and docking using an adaptive artificial potential field,” IEEE 
Transactions on Control Systems Technology, pp. 1-12, 2017, (Under Review) 

R. Zappulla II, H. Park, C. Zagaris, J. Virgili-Llop, and M. Romano, “A metric for the 
equitable comparison of guidance and control methods,” Journal of Guidance, Control, 
and Dynamics, 2017, (Under Review) 

J. Virgili-Llop, J. V. Drew, R. Zappulla II, and M. Romano, “Autonomous capture of a 
resident space object by a spacecraft with a robotic manipulator: Analysis, simulation 
and experiments,” Aerospace Science and Technology, 2017, (Under Review) 


C.4 Conference Proceedings 

R. Zappulla II, H. Park, J. Virgili-Llop, and M. Romano, “Experiments on autonomous 
spacecraft rendezvous and docking using an adaptive artificial potential field approach,” 
in 26th A AS/A!AA SpaceFlight Mechanics Meeting. San Diego, CA: Univelt Inc., 2016, 
paper AAS 16-459, pp. 4461-4478. Available: http://hdl.handle.net/10945/50864 

[Best Paper Award] J. Virgili-Llop, C. Zagaris, H. Park, R. Zappulla II, and M. Ro¬ 
mano, “Experimental evaluation of model predictive control and inverse dynamics con¬ 
trol for spacecraft proximity and docking maneuvers,” 6th International Conference on 
Astrodynamics Tools and Techniques, Mar.14-17 2016, Darmstad, Germany. Available: 
http://hdl.handle.net/10945/50868 

R. Zappulla II, J. Virgili-Llop, H. Park, C. Zagaris, A. Sharp, and M. Romano, “Floating 
spacecraft simulator test bed for the experimental testing of autonomous guidance, 
navigation, & control of spacecraft proximity maneuvers and operations,” AIAA/AAS 
Astrodynamics Specialist Conference, AIAA SPACE Forum, Sep.13-16 2016, AIAA 
Paper 2016-5268. Available: http://dx.doi.Org/10.2514/6.2016-5268 

H. Park, C. Zagaris, J. Virgili-Llop, R. Zappulla II, I. Kolmanovsky, and M. Romano, 
“Analysis and experimentation of model predictive control for spacecraft operations with 
multiple obstacle avoidance,” AIA/VAAS Astrodynamics Specialist Conference, AIAA 
SPACE Forum, Sep.13-16 2016, AIAA 2016-5273. Available: http://dx.d 0 i. 0 rg/l 0.2514/ 
6.2016-5273 

J. Virgili-Llop, J. V. Drew, R. Zappulla II, and M. Romano, “Autonomous capture of a res¬ 
ident space object by a spacecraft with a robotic manipulator: Analysis, simulation and 
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Sep.13-16 2016. Available: http://dx.doi.Org/10.2514/6.2016-5269 

R. Zappulla II, J. Virgili-Llop, and M. Romano, “Near-optimal real-time spacecraft guid¬ 
ance and control using harmonic potential functions and modified RRT*,” 27th AAS/AIAA 
Spaceflight Mechanics Meeting, Feb. 2017, AAS Paper 17-420. Available: http://hdl. 
handle.net/10945/51981 

R. Zappulla II and M. Romano, “A systematic approach to determining the minimum 
sampling rate for real-time spacecraft control,” 27th AAS/AIAA Spaceflight Mechanics 
Meeting, Feb. 2017, AAS Paper 17-424. Available: http://hdl.handle.net/10945/51982 

FI. Park, R. Zappulla II, C. Zagaris, J. Virgili-Llop, and M. Romano, “Nonlinear model 
predictive control for spacecraft rendezvous and docking with a rotating target,” 27th 
AAS/AIAA Spaceflight Mechanics Meeting, Feb. 2017, AAS Paper 17-496. Available: 
http://hdl.handle.net/10945/51981 

J. Virgili-Llop, C. Zagaris, R. Zappulla II, A. Bradstreet, and M. Romano, “Convex op¬ 
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